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CHAPTER  1 
INTRODUCTION 


i.l  Introductory  Comments 

Analysis  of  machine  and  mechanism  problems  prior  to  i960 
relied  upon  various  (and  often  quite  complex)  graphical  techniques. 
These  techniques  were  based  on  geometrical  interpretation  of  con¬ 
straints  and  a  large  body  of  literature  evolved  which  classified 
constraints  and  linkages  [l,  2] .  In  the  1950'"S  digital  and  analog 
computers  began  to  appear,  making  it  possible  to  formulate  and 
repetitively  solve  moderately  sized  systems  of  equations  [5,  . 

This  prompted  the  development  of  mathematical  software  support 
packages  and  the  possibility  of  mathematical  analysis  of  machines 
and  vehicular  systems  became  a  reality.  The  early  programs,  written 
specifically  for  given  applications,  used  geometrical  interpretation 
of  constraints  to  define  relative  position  coordinates ,  thus  helping 
to  reduce  problem  size.  These  programs  were  generally  limited  in 
scope,  requiring  extensive  modification  for  different  applications 
[5,  6]. 

As  computer  speed  and  capacity  increased,  larger  programs 
were  developed,  along  with  data  handling  routines.  It  was  recognized 
that  mechanical  systems  are  composed  of  a  number  of  "standard" 
elements  and  that  they  can  be  combined  as  building  blocks  to  define 
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large  classes  of  mechanisms.  Again,  geometrical  interpretation  of 
constraints  allowed  for  the  reduction  in  the  number  of  independent 
variables,  so  that  considerably  larger,  more  sophisticated  models 
could  be  established.  However,  the  methods  developed  required  top¬ 
ological  analysis  associated  with  identification  of  independent  con¬ 
straint  loops,  which  added  to  program  complexity. 

Several  general  purpose  computer  programs  were  developed  along 
these  lines  in  the  late  60' s  and  early  70’s  [7-9].  These  programs 
are  satisfactory  for  many  mechanisms  applications,  but  incorporation 
of  user  supplied  constraint  and  forcing  functions  is  rather  difficult. 

An  alternate  method  of  formulating  system  constraints  and 
equations  of  motion,  in  terms  of  a  maximal  set  of  coordinates,  by¬ 
passes  topological  analysis  and  provides  for  convenient  user  supplied 
constraints  and  forcing  functions  [10],  This  leads  to  a  more 
general  computer  program,  with  practically  no  limitation  on  mechanism 
type.  The  penalty,  however,  is  a  larger  system  of  equations  to  be 
solved,  which  may  require  greater  computer  capacity  or  time. 

1.2  Program  Efficiency  Versus  Program  Generality 

The  first  general  purpose  computer  codes  capable  of  performing 
dynamic  analysis  of  large  scale,  nonlinear,  constrained,  rigid  body 
mechanical  systems  appeared  around  1970.  The  earliest  programs 
having  major  impact  in  this  area  were  DYMAC  (1970)  [7],  DAMN-DRAM 
(1971)  [8],  IMP  (1972)  [9],  and  ADAMS -3D  (1973)  [10].  Numerous 
programs  have  since  been  developed  but  they  generally  represent 
variations  on  the  above  algorithms  [ll-l4].  These  programs  are 
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similar  in  some  respects  and  are  quite  different  in  others.  Table 
1.1  lists  some  of  the  major  features  of  the  above  four  programs. 

In  the  design  of  algorithms  it  is  intended  that  program 
efficiency  and  generality  be  maximized.  Efficiency  is  easy  to 
measure.  Generality,  however,  is  highly  dependent  upon  the  needs  of 
the  user  and  a  program  that  is  ideally  suited  for  certain  applications 
may  be  completely  inadequate  for  others.  Therefore  it  is  difficult 
to  assign  ratings  tel  generality  of  a  given  code. 

Many  codes  are  developed  using  linkages  as  subelements  of 
mechanical  systems.  Since  linkages  generally  form  loops  (closed  or 
open)  a  minimum  number  of  linkage  relative  position  coordinates 
(designated  as  Lagrangian  coordinates  by  Paul  [15])  are  defined, 
allowing  loop  constraint  equations  to  be  formulated.  Complex  systems 
may  have  many  loops  (not  all  independent)  so  graph  theory  or  user 
preprocessing  is  required  to  avoid  unnecessary  formulation  of  redun¬ 
dant  loop  constraint  equations  [16,  17].  Loop  constraint 
formulation  and  topological  preprocessing  requirements  add  an  order 
of  magnitude  of  complexity  to  codes  when  the  user  desires  to  incor¬ 
porate  constraint  relations  that  are  not  provided  in  the  standard 
formulation.  A  definition  of  program  generality  might  then  be  a 
measure  of  how  easily  one  supplements  the  standard  code  to  provide 
for  unanticipated  constraint  or  forcing  elements. 

Constraint  loop  formulations,  with  topological  preprocessing, 
are  desirable  from  a  numerical  standpoint  since  they  allow  formula¬ 
tion  of  a  problem  using  a  minimum  number  of  state  variables.  They 
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lead,  to  a  corresponding  minimum  number  of  highly  coupled  constraint 
equations.  The  constraint  Jacobian  matrix  resulting  from  lineari¬ 
zation  of  constraint  equations  will  generally  be  small,  even  for 
relatively  large  problems.  Well  established  f'lll  matrix  manipulation 
algorithms  can  be  efficiently  applied.  The  small  number  of  state 
variables  thus  yields  a  small  number  of  differential  equations  of 
motion,  which  is  a  definite  advantage  for  numerical  integration 
algorithms  [18]. 

The  programs  DYMAC,  DRAM,  and  IMP  employ  loop  methods,  as 
described  above.  A  different  approach,  however,  is  taken  in  ADAMS 
[10,  19J  whereby  constraint  equations  are  formulated  between 
bodies  connected  by  given  joints.  The  concept  of  lin/, '  is  not  used 
in  this  development.  Elements  of  the  system  are  treated  as  rigid 
bodies  each  with  a  Cartesian  coordinate  system  attached.  One  does 
not  think  in  terms  of  element  type,  i.e.,  slider,  crank,  rocker, 
etc.;  but  in  terms  of  joint  type,  i.e.,  revolute,  translational, 
cylindrical,  etc.  Since  each  system  element  has  an  identical  coor¬ 
dinate  system,  arbitrary  user -supplied  constraints  and  generalized 
forcing  functions  can  be  formulated  without  regard  to  element  type. 
Thus,  equation  formulation  and  computer  programming  are  simplified. 

A  major  disadvantage  of  this  method  arises  from  initially 
assigning  a  maximal  set  of  generalized  coordinates  (degrees  of 
freedom)  to  the  system.  Nonlinear  algebraic  constraint  equations 
are  then  imposed  to  remove  system  degrees  of  freedom.  This  results 
in  a  maximal  number  of  differential  and  algebraic  equations,  which 
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traditionally  have  been  solved  iteratively  and  simultaneously  by 
stiffly  stable  implicit  numerical  integration  algorithms  [20]. 

Solving  large  systems  of  equations  iteratively,  even  when  advantage 
is  taken  of  matrix  sparsity,  is  time  consuming.  In  addition,  the 
corrector  equations  thus  obtained  contain  variables  related  to 
accelerations,  on  which  error  control  cannot  be  maintained.  Poor 
prediction  of  these  variables  results  in  an  excessive  number  of 
corrector  iterations  (up  to  five  or  more),  forcing  even  smaller  time 
steps  than  would  otherwise  be  necessary. 

In  summary,  the  loop  method  of  constraint  formulation  gains 
program  efficiency  at  the  expe  ;e  of  program  generality,  whereas  the 
joint  method  of  constraint  formulation  in  ADAMS  gains  program  gener¬ 
ality  at  the  expense  of  program  efficiency.  The  objective  of  the 
research  reported  in  this  report  is  development  of  a  numerical 
analysis  method  that  contains  the  strong  points  of  both  these  methods. 

1.3  Obtaining  Efficiency  Without 
Sacrificing  Generality 

The  main  objective  of  this  research  is  to  develop  numerical 
and  analytical  methods  to  improve  the  efficiency  of  general  purpose 
analysis  programs  that  determine  dynamic  response  of  large  scale, 
constrained,  rigid-body  mechanical  systems.  A  planar  dynamic  analysis 
code  [10,  19,  21],  similar  to  ADAMS,  formulates  nonlinear  holo- 
ncraic  constraint  equations  and  differential  equations  of  motion, 
written  in  terms  of  a  maximal  set  of  Cartesian  generalized  coordi¬ 
nates,  to  facilitate  the  general  formulation  of  constraints  and 
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forcing  functions.  A  maximal  set  of  coupled  algebraic  and  differ¬ 
ential  equations  are  then  solved  iteratively  and  simultaneously  by 
an  implicit  numerical  integration  algorithm  [ 20] .  As  noted  earlier 
the  method  is  inefficient  (even  though  sparse  matrix  manipulation 
algorithms  are  employed) ,  because  an  excessive  number  of  equations 
are  involved  in  the  iteration  process  within  the  numerical  integra¬ 
tion  algorithm. 

A  modified  approach  developed  in  this  report  maintains  program 
generality  by  keeping  the  above  method  of  constraint  and  equation  of 
motion  formulation.  An  intermediate  numerical  processing  step  is 
introduced  that  effectively  eliminates  the  constraint  equations  and 
dependent  equations  of  motion,  before  each  numerical  integration 
step.  Iteration  is  then  limited  to  solving  only  the  set  of  constraint 
equations.  The  time  saved  in  the  iteration  step  and  in  the  numerical 
integration  step  far  outweighs  the  overhead  introduced  by  the  inter¬ 
mediate  numerical  processing  step.  This  method  has  demonstrated  an 
order  of  magnitude  increase  in  speed,  when  highly  constrained  systems 
are  simulated. 

A  Gaussian  elimination  algorithm  with  full  pivoting  [22] 
decomposes  the  constraint  Jacobian  matrix  and  identifies  dependent 
and  independent  generalized  coordinates.  For  small  systems,  the 
algorithm  constructs  an  influence  coefficient  matrix  relating  varia¬ 
tions  in  dependent  and  independent  variables.  For  large  systems,  in 
which  matrix  sparsity  becomes  a  significant  factor,  tlv*  algorithm 
provides  information  necessary  to  set  up  a  modified  sparse  matrix 
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that  relates  variations  in  dependent  and  independent  variables.  In 
either  case,  this  information  is  employed  to  numerically  construct  a 
reduced  system  of  differential  equations  of  motion,  whose  solution 
yields  the  total  system  dynamic  response. 

A  pieced  interval  analysis  concept  is  also  developed  to 
further  the  enhancement  of  program  efficiency.  When  rapidly  changing 
or  "intermittent"  events  occur  within  a  simulation  interval,  it  may 
be  desirable  or  necessary  to  break  the  time  interval  into  subintervals. 
Different  equations  or  analysis  techniques  may  then  be  incorporated 
at  these  subintervals.  For  example,  numerical  integration  algorithms 
employing  high  order  polynctnial  predictors  become  inefficient  or  fail 
near  abrupt  or  discontinuous  variations  of  state  variables.  When 
event  times  are  known  in  advance,  time  stepsize  control  can  be  placed 
on  the  integration  algorithm  before  encounter,  to  avoid  lengthy 
search  for  correct  reduced  time  step  sizes.  Momentum  balance  is 
performed  to  determine  jump  discontinuity  in  velocity  if  bodies 
impact  or  impulsive  loading  occurs. 

Event  occurrences,  such  as  impact  between  bodies,  are 
generally  functions  of  system  state  or  time.  Therefore  an  event 
predictor  is  incorporated  into  the  integration  algorithm.  Logical 
functions  of  state  and  time  are  formulated  such  that,  as  each  given 
function  passes  through  zero,  it  defines  a  logical  time  at  which 
other  actions  may  be  taken.  Logical  functions  are  extrapolated  ahead 
before  the  system  state  is  advanced,  so  that  significant  events  may 


9 


be  detected  before  being  encountered  by  the  system.  The  logical  time 
corresponding  to  a  zero  of  a  logical  function  is  determined  and  the 
system  solution  is  forced  precisely  at  this  time. 

l±k _ Scone  of  the  Report 

Both  planar  and  three  dimensional  rigid-body  dynamic  analysis 
programs,  employing  Cartesian  coordinates  and  a  simple  constraint 
formulation,  have  been  developed  [10,  19,  21,  23,  24],  In  addition, 
flexible  beam  linkage  elements  [25]  have  been  incorporated  into 
the  two  dimensional  program  to  extend  its  range  of  applications. 

The  analysis  methods  developed  in  this  research  apply  equally  well 
to  each  of  the  above  programs.  In  the  interest  of  economy  only  the 
planar  rigid-body  analysis  method  is  developed  here.  This  is  suf¬ 
ficient  to  illustrate  the  analysis  procedure,  without  introducing 
unnecessary  complexity. 

In  Chapter  2,  the  planar  rigid-body  dynamic  analysis  algorithm 
is  developed.  An  implicit,  stiffly  stable  numerical  integration 
algorithm  is  employed  to  simultaneously  and  iteratively  solve 
algebraic  constraint  equations  and  differential  equations  of  motion. 
In  addition,  it  is  shown  that  this  method  of  solution  results  in  a 
sparse  Jacobian  matrix  of  maximal  dimension,  and  that  sparse  matrix 
manipulation  algorithms  must  be  employed  to  gain  efficiency,  even  for 
small  problems.  Further,  it  is  shown  that  the  method  requires  an 
excessive  number  of  corrector  iterations,  employing  the  above  matrix 
at  each  time  step. 
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In  Chapter  3  the  "inverse  dynamics  problem"  is  presented 
[26].  That  is,  assuming  that  motion  of  certain  generalized  coor¬ 
dinates  (equal  in  number  to  system  degree  of  freedom)  is  specified, 
the  entire  system  state,  including  all  external,  inertial,  and  con¬ 
straint  forces  can  be  determined.  This  is  done  by  an  iterative 
algebraic  procedure,  requiring  no  numerical  integration.  Methods 
for  reducing  problem  size  by  numerical  elimination  of  constraint 
equations  and  dependent  variables  are  more  easily  demonstrated  for 
the  inverse  dynamics  problem.  However,  these  methods  are  equally 
applicable  to  the  dynamics  problem  of  Chapter  4. 

It  is  shown  that  the  Jacobian  matrix  formed  by  linearizing  the 
constraint  equations  expresses  relations  between  variations  in  the 
generalized  coordinates.  The  above  matrix  also  expresses  similar 
relations  between  higher  order  derivatives.  A  Gaussian  elimination 
procedure  is  employed  to  identify  a  submatrix  of  maximal  rank  and  it 
is  shown,  by  the  implicit  function  theorem,  that  generalized  coor¬ 
dinates  corresponding  to  certain  columns  of  this  matrix  depend 
entirely  on  the  remaining  coordinates  and  possibly  time,  hence  they 
are  designated  dependent  variables.  These  dependent  variables  may 
be  evaluated  if  the  independent  variables  are  known,  since  this 
submatrix  is  nonsingular.  Likewise,  dependent  velocities  and  accel¬ 
erations  are  evaluated,  when  corresponding  independent  velocities 
and  accelerations  are  known.  Constraint  reaction  forces,  included 
in  the  equations  of  motion,  are  given  by  the  product  of  the  transpose 
of  the  above  matrix  and  a  vector  of  seeder  multipliers.  Eliminating 
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the  scaler  multipliers,  employing  nonsingularity  of  the  above  sub- 
matrix,  yields  a  minimal  set  of  equations  of  motion,  involving  only 
independent  variables  and  excluding  all  constraint  reaction  forces. 

A  significant  number  of  intermediate  matrix  products  are 
involved  in  the  above  procedure.  For  large  systems,  efficiency  would 
be  lost  when  carrying  out  the  indicated  sparse  matrix  products. 
Therefore,  an  alternate  sparse  matrix  method  is  developed  that 
requires  no  intermediate  matrix  products,  yet  effectively  arrives  at 
the  same  minimal  set  of  equations  of  motion. 

An  additional  problem  arises  when  solving  nonlinear  equations 
by  an  iterative  procedure,  such  as  Newton's  method.  If  the  initial 
estimate  of  roots  of  a  system  of  nonlinear  equations  is  poor,  con¬ 
vergence  will  be  slow  or  may  fail.  While  the  methods  developed  in 
Chapter  3  minimize  the  effect  of  errors  in  dependent  variables  (on 
convergence  rate),  better  estimates  will  result  in  fewer  or  quite 
often  no  iterations .  When  advancing  a  constrained  system  from  one 
>  position  to  another,  if  dependent  velocities  and  accelerations  are 

determined,  they  may  be  used  to  arrive  at  better  extrapolated  values 
for  dependent  position  variables.  Alternately,  a  variable  order 
polynomial  extrapolator  fitted  to  dependent  variable  history  will 
yield  even  better  results.  Newton's  backward  divided  difference 
|  formula  is  very  convenient  and  efficient  for  this  application  [27]. 

The  numerical  procedures  of  Chapter  3  are  extended  to  the 
direct  dynamics  problem  in  Chapter  4  [26].  Here  it  is  assumed 
that  one  or  more  of  the  system  generalized  coordinates  are  not  given 
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as  functions  of  time.  Therefore,  the  equations  of  motion  must  he 
solved  by  numerical  integration  to  determine  system  response  to 
forces.  The  minimal  system  of  differential  equations  of  motion  may 
then  be  conveniently  solved  by  explicit  predictor/implicit  corrector 
multistep  methods,  such  as  the  Adams  Bashf orth/Adams  Moulton  method 
[28]. 

Additional  problems  arise  as  a  result  of  identifying  optimal 
sets  of  independent  variables.  The  sets  generally  change  as  system 
configuration  changes.  Therefore,  the  minimal  set  of  differential 
equations  also  changes.  To  allow  for  the  possibility  that  any  vari¬ 
able  may  become  independent,  polynomial  extrapolators  are  maintained 
for  all  dependent  variables.  A  given  polynomial  extrapolator  auto¬ 
matically  charges  to  a  predictor  if  the  corresponding  variable 
becomes  independent,  and  a  predictor  reverts  back  to  an  extrapolator 
when  a  variable  becomes  dependent.  This  procedure  provides  accurate 
estimates  for  the  Newton  iteration  described  earlier  for  solution  of 
constraint  equations  and  avoids  interruption  of  the  numerical  inte¬ 
gration  procedure  caused  by  lack  of  variable  history. 

The  concept  of  pieced  interval  analysis  is  introduced  in 
Chapter  5.  In  the  interest  of  efficiency  and  simplicity  a  given 
simulation  may  be  divided  into  two  or  more  time  intervals  for  which 
different  governing  equations  or  analysis  procedures  may  apply. 
Boundaries  of  such  intervals  may  correspond  to  instances  where 
violent  actions  such  as  member  impact,  impulsive  loading,  mass  cap¬ 
ture,  mass  release,  member  breaking,  etc.  occurs.  Such  events, 
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when  treated  in  a  continuous  manner,  create  serious  problems  for  high 
order  predictor/corrector  integration  algorithms,  when  their  presence 
is  not  known  in  advance  [29J.  The  ability  to  force  small  step  sizes 
and/or  relax  error  requirements  just  prior  to  event  encounter  is  an 
effective  means  of  improving  program  efficiency.  If  member  distortion 
is  significant,  a  temporary  flexible  model  may  be  implemented.  Times 
at  which  such  events  occur  are  usually  not  known  in  advance,  since 
they  depend  on  system  state.  Therefore,  logical  functions  of  system 
state  and  time,  supplied  by  the  user,  are  employed  by  the  integration 
algorithm  to  predict  event  occurrences  before  encounter.  These 
logical  event  occurrences  define  logical  times  that  mark  the  interval 
noundaries . 

A  logical  events  monitor  is  developed  in  Chapter  5*  In  addi¬ 
tion,  momentum  balance  equations  for  general  planar  constrained  rigid 
body  systems  undergoing  impulsive  loading  or  impact  are  derived. 
Incorporation  of  flexible  degrees  of  freedom  into  the  model  has 
been  dealt  with  [25],  and  will  not  be  included  in  this  report.  It 
is  interesting  to  note  that  if  a  flexible  model  is  required  on  an 
interval  for  which  small  angular  rotations  occur,  the  problem  is 
essentially  linear  and  efficient  modal  analysis  techniques  may  be 
employed. 

Numerical  examples  are  presented  in  Chapter  6  to  demonstrate 
efficiency  and  generality  of  the  method.  However,  a  detailed  descrip¬ 
tion  of  the  operational  and  data  requirements  of  the  computer  program 
is  not  presented  in  this  report.  This  information  will  be  available 
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in  a  theoretical  and  user's  manual  that  will  allow  periodic  update, 
since  much  of  the  program  is  still  in  a  developmental  stage.  Examples 
are  presented  in  enough  detail  to  demonstrate  the  contributions  of 
this  research.  In  addition,  they  suggest  applications  not  found  in 
general  purpose  mechanical  systems  analysis  programs. 

The  first  example  treated  investigates  dynamic  pitch  response 
of  an  articulated  vehicle,  consisting  of  two  tracked-vehicles  that 
are  coupled  together  by  an  articulation  joint  and  a  hydraulic 
actuator.  A  representative  model  of  the  mechanical  vehicular  system 
is  developed,  primarily  from  the  standard  program  element  formula¬ 
tion.  Additional  user  supplied  generalized  forcing  functions  and 
constraints  are  formulated  to  incorporate  major  contributions  to 
system  dynamic  response.  An  electro -mechanical -hydraulic  control 
system,  with  position  and  force  feedback,  is  formulated  and  coupled 
to  the  mechanical  system  through  extensions  of  a  standard  element  in 
the  program.  This  system  controls  relative  vehicular  pitch  attitude 
by  monitoring  hydraulic  actuator  extension  rate  and  system  pressure. 
The  control  system  differential  equations  and  mechanical  system 
differential  equations  of  motion  are  solved  simultaneously  by  the 
same  integration  algorithm.  Thus,  no  iterative  procedures  or 
approximations  are  required  to  obtain  the  coupled  system  dynamic 
response. 

/ 

A  second  example  illustrates  the  pieced  interval  analysis 


capability  of  the  program.  Recoiling  dynamics  of  a  weapon  system, 
subject  to  discontinuous  and  irrqpulsive  forces,  impacts,  and  mass 
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capture  and  release  is  investigated.  In  addition,  the  system  con¬ 
tains  a  cammed  loading  mechanism  that  affects  barrel  dynamics. 

Logical  functions  are  developed  for  this  intermittent  motion  problem. 
It  is  shown  how  various  programmed  actions  are  taken,  in  response  to 
each  active  logical  function.  For  this  problem,  three  actions  in 
order  of  complexity  are: 

1.  Remove  or  add  forces  to  the  system.  Restart  the  integra¬ 
tion  process. 

2.  Apply  impulsive  loads  to  the  system.  Perform  momentum 
balance  to  determine  jump  discontinuities  in  velocity. 
Restart  the  integration  process. 

3.  Impact,  mass  capture,  or  mass  release  occurs.  Introduce 
restitution  equations  and  perform  momentum  balance  to 
determine  jump  discontinuities  in  velocity.  Impose  or 
remove  constraints  for  mass  capture  or  release,  respec¬ 
tively.  Restart  the  integration  process. 
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CHAPTER  2 

THE  PLANAR  RIGID  BODY  MODELING  METHOD 
FOR  CONSTRAINED  SYSTEMS 

2.1  Introduction 

The  term  "Rigid  Body  Systems"  distinguishes  the  intent  of  the 
analysis  method  and  program  developed  in  this  chapter.  Many  of  the 
analysis  programs  existing  today  are  based  on  the  traditional  notion 
of  "machines  and  mechanisms,"  in  which  the  ideas  of  "links"  and 
"closed  loop  analysis"  are  embedded  in  the  "geometry  of  mechanisms". 
This  limitation  complicates  an  orderly  approach  to  general  constrained 
dynamic  system  analysis . 

Lagrange  may  have  been  the  first  to  recognize  the  limitation 
of  the  geometrical  approach,  as  evidenced  by  the  statement  found  in 
the  preface  of  his  Mecanique  Analytique  [50]:  "No  diagrams  will  be 
found  in  this  work.  The  methods  that  I  explain  in  it  require  neither 
constructions  nor  geometrical  or  mechanical  arguments,  but  only  the 
algebraic  operations  inherent  to  a  regular  uniform  process.  Those 
who  love  Analysis  will,  with  joy,  see  mechanics  become  a  new  branch 
of  it  and  will  be  grateful  to  me  for  thus  having  extended  its  field." 

The  concept  of  constraints  in  physical  systems  is  not  unlike 
the  concept  of  circuits  in  electrical  systems.  Kirchhoff  stated 
two  basic  laws  that:  (l  -  voltage  law)  "The  voltages  with  their 
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proper  signs,  taken  completely  around  a  mesh,  add  up  to  zero",  and 
(2  -  current  law)  "The  currents  at  a  node,  with  due  regard  to 
direction,  add  up  to  zero."  From  these  laws  and  Ohm's  law,  two  basic 
methods  of  circuit  solution  have  developed  as  follows: 

(a)  the  mesh  method  of  solution,  characterized  by  the  sum¬ 
mation  of  voltages  around  meshes,  is  an  application  of 
the  voltage  law. 

(b)  the  nodal  method,  characterized  by  the  summation  of 
currents  at  nodes,  is  an  application  of  the  current  law. 

Loop-closure  methods  of  mechanism  analysis  are  analogous  to 
the  mesh  method,  (a).  An  alternate  and  more  direct  approach,  first 
developed  on  a  large  scale  by  Orlandea  [lO]  for  general  three 
dimensional  analysis,  is  based  on  the  nodal  method,  (b).  He  called 
this  method  "the  Node  -  Analogous  Formulation"  for  mechanical  systems 
derived  by  considering  that  D'Alembert's  Principle  for  forces  and 
Kirchhoff's  law  for  currents  are  analogous.  He  further  identified  a 
finite  number  of  standard  components  of  mechanical  systems  (linkages) 
While  his  method  easily  allows  all  of  the  constraints  associated  with 
linkages,  there  need  be  no  connotation  of  linkages  in  its  development 

The  analysis  method  presented  here  is  developed  according  to 
the  nodal  representation,  as  introduced  by  Orlandea.  Since  this 
method  is  nearly  devoid  of  geometrical  and  topological  concepts,  it 
is  adapted  to  a  broad  class  of  mechanical  systems.  These  systems 
may  range  from  the  most  complex  linkage  mechanisms  to  articulated 
tracked  vehicles  with  electro-hydraulic  control  systems. 


18 


In  this  chapter,  a  general  system  of  constrained  equations  of 
rigid  body  planar  motion  is  formulated.  An  implicit  numerical  inte¬ 
gration  method  for  automatically  solving  the  system  equations  and 
sparse  matrix  techniques  are  also  discussed.  Generalized  coordinates 
are  defined  to  locate  each  individual  rigid  member  of  the  system  and 
to  express  kinetic  energy  for  each  member.  Constraints  between 
elements  are  taken  as  friction -free  standard  constraints,  with  pro¬ 
visions  for  additional,  non-standard  constraint  formulation,  as 
needed.  In  addition  to  standard  constraints,  springs,  dampers,  and 
actuators  connecting  any  pair  of  points  on  different  bodies  of  the 
system  are  included  in  the  model.  These  standard  force  elements, 
together  with  allowance  for  arbitrary  non-standard  forcing  functions, 
make  the  formulation  quite  general. 

Implicit  numerical  integration  and  sparse  matrix  algorithms 
C31,  32]  were  initially  used  in  numerical  integration  of  the 
equations  of  motion.  After  considerable  numerical  experimentation, 
it  became  apparent  that  implicit  simultaneous  solution  of  algebraic 
and  differential  equations  create  artificially  stiff  [20]  problems 
and  the  integration  of  large  systems  of  equations  increases  the 
spectrum  of  frequencies  with  which  integration  algorithms  must  cope. 
This  prompted  a  search  for  new  methods  of  solving  the  system  equa¬ 
tions  of  motion,  which  forms  the  basis  of  the  method  presented  in 
the  following  chapters.  The  implicit  method  is  presented  in  this 
chapter,  for  comparison  purposes  only.  The  basic  constraint 
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formulation  and  equations  of  motion  presented  in  this  chapter  are 
used  as  the  basis  for  developments  in  the  following  chapters. 

2.2  Constrained  Equations  of  Planar  Motion 

(a)  Generalized  Coordinates  and  Energy  Equations 

In  order  to  determine  the  configuration  or  state  of  a  planar 

mechanical  system,  it  is  first  necessary  to  define  generalized 

coordinates  that  specify  the  location  of  each  body  in  the  system.  As 

shown  in -Fig.  2.1,  let  the  x-y  coordinate  system  be  a  fixed  inertial 

reference  frame.  Define  a  body -fixed  ^  -  T]^  coordinate  system 

embedded  in  a  typical  body  i.  The  location  of  body  i  is  specified 

by  the  global  coordinates  (x^,  y^)  or  the  vector  of  the  origin  of 

its  reference  frame  and  the  angle  cp  of  rotation  of  the  body  fixed 

5^ -axis,  relative  to  the  global  x-axis. 

The  center  of  mass  of  body  i  is  located  by  a  vector  r  in 

i 

the  body  fixed  coordinate  system,  with  -  T!^  components  and 

i 

T|  .  In  terms  of  the  ger  ''ralized  coordinates  x.  ,  y,  ,  and  and  the 
nu  i*  Ji'  Yi 

parameters  £  and  T]  ,  the  center  of  mass  is  located  globally  by 
i  i 


x  =  x.  +  t  cos  ra.  -  T  sm  <p. 

m.  i  sm.  ri  m.  Yi 

i  i  i 

y  =  y.  +  P  sin  m.  +  T1  cos  cp. 

m.  i  Yi  m.  Yi 

i  i  i 


(2.1) 


• ... *  ■ 
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Thus,  the  kinetic  energy  of  body  i  is 


(2.2) 


where  m  is  the  mass  of  body  i  and  is  its  centroidal  polar  moment 
of  inertia. 


(b)  Equations  of  Constraint 

Figure  2.1  further  depicts  an  adjacent  body  j,  with  body- 

fixed  coordinate  system  located  by  the  vector  $ . .  Let  arbitrary 

points  p. .  on  body  i  and  p..  on  body  j  be  located  by  vectors  . 

ij 

and  r..,  specified  in  the  body  fixed  coordinate  systems  by  coordi- 

nates  ,  T^,  5^,  and  T|^.  These  points  are  in  turn  connected  by 

a  vector  r  , 

P 


r  =  +  r.  .  -  1$  -  r . . 

P  i  ij  j 


(2.3) 


The  vector  condition  for  a  rotational  joint  between  bodies  i 

and  j,  at  points  p. .  and  p..,  is  simply  r  =  3,  yielding  the  follow- 

J  1  P 

ing  pair  of  scaler  constraint  equations: 


x.  *  5y  cos  -  q.  sin  cos  Vj 


+  hi  Sl"  ”3  ‘  ° 


22 


yi  +  ?ij  Sln  +  \j  C°S  'Pi  '  yj  "  ?ji  Sin  'Pj 

-  T| . .  cos  cp  .  =  0  (2.4) 

If  r^  in  Fig.  2.1  is  taken  as  a  nonzero  vector  of  fixed  length 

r  ,  a  "massless  link"  of  length  r  connects  points  p. .  and  p .. .  A 
p  P  ij  Ji 

single  scalar  equation  for  this  constraint  can  then  be  written  as 


V  sln  h)2  * 


(yi  + 


«ij  sln  fl 


COS 


yi  '  !di  Sln  tl  "  V  °OS 


r2  -  0 
P 


(2.5) 


For  a  translational  joint  shown  in  Fig.  2.2, 

and  p  on  body  i,  and  points  p  and  p  on  body  j 
J-*-  J  -L  J 

parallel  to  the  path  of  relative  motion  between  the 

that  the  specified  body  fixed  vectors  and  are 

nitude.  Since  1?,  and  ^  are  parallel,  x  ^  =  "5, 

J  ^  J 

component  yielding  the  scalar  equation 


let  points  p  ^ 
lie  on  some  line 
two  bodies ,  such 
of  nonzero  mag- 
with  zero  z 


[(?i2  ‘  cos  <Pi  '  (\2  "  11il)  Sin  *PiJ 

[(5j2  -  5jl>  -  V  '°s  ^ 
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'  5jl)  COS  'Pj  ‘  (T1J2  "  Sin  'Pj ] 

*  ^ii.)  ^Pj_  +  (^1^2  "  3  =  0  (2.6) 

Likewise ,  r..  and  "3.  are  parallel,  so  r . .  x  yields  a  zero  z 

J  Ji  j 

component  and  the  second  scalar  equation 


[xi  +  ?il  COS  'f’i  "  \l  sin  Vi  -  xj  -  5j!  cos  <Pj  +  ■njl  sin  cp  J 
[(?J2  *  5jl)  Sln  'Pj  +  (Tlj2  "  T1tjl)  COS  'Pj] 

-  [y±  +  ?il  sin  <pi+n.1  cos  cp.  -  y^  -  f.±  sin  cp^  -  T]^  cos  <p.] 
[(?j2  ‘  5ji)  cos  -  (Hj2  -  H^)  sin  q>  ]  =  0  (2.7) 


The  parameters  (g^,  7lil)  and  (|i2>  l^)  locate  points  p  and  pig 
in  body  i  coordinate  system,  and  (gjr  T)^)  and  (5^,  7)  )  locate 

points  p^  and  p^2  in  body  .j  coordinate  system. 


Other  constraints  may  be  formulated  by  a  similar  process, 
i  T 

Denote  by  q  =  [x  ,  y  ,  cp  ]  the  vector  of  generalized  coordinates 
1  1  T1  5T  T 

of  body  i  and  by  q  =  [q  ,  q"  ,  • • . ,  q°  the  composite  vector  of 
all  system  generalized  coordinates.  In  this  notation,  the  holonomic 
constraints  of  Equations  2.4  to  2.7  and  other  (perhaps  time  dependent) 
holonomic  constraints  can  be  written  in  vector  function  form  as 
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«(q,t)  =  0 


(2.8) 


T 

where  $(q,t)  =  [$  (q,t),  ...,  §m(q,t)j  and  the  functions 

i  =  1,  . . . ,  m,  arc  assumed  to  be  independent,  i.e.,  the  Jacobian 
9?/dq  =  [S5./3q.]  has  full  row  rank. 

(c)  Equations  of  Springs,  Dampers 
and  Actuators 

Internal  forces  due  to  other  types  of  elements  acting  between 
bodies  may  be  obtained  by  a  process  similar  to  the  constraint 
equation  development.  For  example,  since  springs,  dampers,  and 
actuators  shown  in  Fig.  2.3  generally  appear  together,  they  are 
incorporated  into  a  single  set  of  equations.  The  equation  for 
spring -damper -actuator  force  is 


?.  .  =  k.  .(£.  .  -  Zn  )  +  c.  .  i.  .  +  F„ 

«  L  °ij  hd  sy 


where  .  is  the  resultant  force  vector  F  l  +  F,,  j  in  the 
10  x.  yij 

-f  ° 

element,  R  is  the  vector  Z.  ■  cos  a  l  +  Z.  .  sin  a  o  between  points 
s .  .  10  10 

J 

s.  .  and  s...  k.  .  is  an  elastic  spring  coefficient  that  may  depend  on 
10  01  10 

generalized  coordinates  and  time,  c„  is  a  damping  coefficient  that 

may  depend  on  generalized  coordinates  and  time,  l  is  the  undeformed 

°ij 

spring  length,  l  is  the  deformed  spring  length,  Z.  .  is  the  time 
i*0  iJ 

derivative  of  Z.  .,  and  F_  is  an  actuator  force  applied  along  the 

^  ij 

element  that  may  depend  upon  generalized  coordinates  and  time.  The 
unit  vectors  1  and  3  are  parallel  to  the  x  and  y  axes,  respectively. 


?6 
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(d)  Conservative  and  Nonconservative 
Generalized  Forces 

Contributions  of  forces  acting  on  body  i  to  the  system  equa¬ 
tions  of  motion  are  determined  by  employing  a  virtual  work  concept. 
The  virtual  work  of  externally  applied  forces  and  spring -damper - 
actuator  forces  acting  on  body  i  is  written  as 

fW1  =  Q.^(q,q,t)  6x  +  Q^(q,q,t)  6y.  +  Q.^(q,q,t)  by.  (2.10) 

j  vr 

x  l  i.  i  T 

The  vector  Q  =  [G  ,  Q  ,  Q  ]  of  generalized  forces  on  body  i  is 
x  y  T  cp  T  T 

12  n  T 

thus  defined  and  Q  =  [0  ,  Q  ,  . . . ,  Q  ]  is  the  vector  of  system 

generalized  forces.  Typical  generalized  forces  are  calculated  to 
illustrate  the  procedure.  Consider  here  the  spring -damper -actuator 
element  of  Fig.  2.3  where  point  s. .  on  body  i  is  located  by  the 
vector  +  r  .  A  virtual  change  in  the  location  of  point  s. .  is 

1  Sij  10 

given  by 


Forming  the  dot  product  of  f.  ,  with  6  ($.  +  r  )  determines  contri- 

ij  i  8^ 


butions  to  the  generalized  force  expressions  for  body  i  as 
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Generalized  force  contributions  are  summed  for  all  externally  and 
internally  applied  forces  acting  on  body  i,  to  arrive  at  Q1. 

(e)  System  Equations  of  Motion 
Virtual  displacements  6q  that  are  consistent  with  constraints 
(i.e.,  with  time  fixed)  satisfy 

5q  6q  =  0  (2.13) 

where  5q  =  d?/dq  =  [SS^/dq^ using  subscript  notation  for  dif¬ 
ferentiation  with  respect  to  a  vector. 

The  variational  form  of  Lagrange's  equations  of  motion,  with 
workless  constraints,  is  [33] 


(2.1*0 
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which  must  hold  for  all  6q  satisfying  Equation  2.13.  By  Parkas' 

Lemma  (3U],  there  exists  a  vector  k  e  Rm  of  multipliers  (called 
Lagrange  multipliers )  such  that 

fr  (t-)T  -  (T  )T  -  Q  -  A  =  o  (2.15) 

dt  q  q  q 

which  with  Equation  2.8  form  the  constrained  equations  of  motion  of 
the  system.  For  planar  systems  treated  here.  Equations  2.1  and  2.2 
yield 


q 


T 


(T1i)T  =  M1  (q1)  q1  -  A1  (q\  q) 


(2.16) 


where 


S 


K1  - 


0 


|^-®1(SBi  »in  9t  *  COB  9^ 


0  l,ln  *1  *  \  COB  V 

°i  -“i(*  \  co*  ’i  +  \  sln  V 

-V-  \  co‘  *  \  ,ln  ’i5  Ji  *  "i  ^ 


(2.17) 


50 


and 


m 


.2 


t.  (f  cos  tp.  —  T)  sin  cp.  )cp. 
1  Tn.  Yi  m.  Yi  Yi 
l  i 


m.  (?  sin  m.  +  T|  cos  cp.  )cp 
l  sm.  Yi  m.  Yi  Y: 


.2 


S 


Using  matrix  notation,  M  =  d±ag(M\  M^,  M11) and 

.  T  2^  t 

A  =  [A  ,  A  ,  A  ]  ,  Equation  2.15  becomes 


M(q)q  =  A(q,q)  +  Q(q,q,t)  +  Sq(q,t)X 


Initial  conditions  for  system  motion  are  given  as 

q(tQ)  =  q° 
q(t0)  =  q° 

where  q°  and  q°  are  consistent  with  constraints,  i.e.,  q 
Equation  2.8  and  q*^  satisfies 


a  q  +  a  = 
q  t 


(2.18) 


(2.19) 


(2.20) 

satisfies 


0 


(2.21) 
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2.3  Integration  of  Constrained  Equation 
of  Motion  by  an  Implicit  Method 

(a)  Numerical  Solution  of  Algebraic  and. 

Differential  Equations 

The  method  of  numerically  integrating  differential  equations 
presented  in  this  and  following  chapters  requires  that  they  be 
reduced  to  first  order  form.  This  is  accomplished  by  introducing 
the  vector 


w  =  q 


into  Equation  2.19  to  obtain 


M(q)w  -  A(q,w)  -  Q(q,w,t)  -  5^(q,t)\  =  0 


with  initial  conditions 
q(tQ)  =  q° 
w(tQ)  =  w°  =  q° 


where 


?q(q°,t0)  w  (tQ)  +  ?t(q°,t0)  =  0 


(2.22) 


(2.23) 


(2.24) 


(2.25) 
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In  order  to  solve  the  differential  equations  of  motion, 
numerical  integration  theory  is  used  to  obtain  a  set  of  approximations 
that  is  suitable  for  digital  computation.  Integration  of  the  equa¬ 
tions  of  motion  is  accomplished  as  a  simultaneous  solution  of 
algebraic  and  differential  equations.  Standard  approaches,  however, 
are  designed  to  solve  systems  of  differential  equations  of  the  form 

y  =  f(y,t)  (2.26) 

where  y  is  an  m-vector  variable  and  f  is  an  m-vector  of  functions. 

A  modified  approach  is  taken  here  that  allows  for  the  simultaneous 
solution  of  mixed  algebraic  and  differential  equations  of  the  form 

g(y,y,t)  =  0  (2 .27) 

where  y  may  not  appear  in  some  of  the  equations.  Before  introducing 
the  method  to  be  used,  it  is  instructive  to  review  the  standard 
approach  applied  to  solving  Equation  2.26. 

The  basic  method  of  constructing  approximate  solutions  is  to 
place  a  grid  of  time  points  t^,  i  =  1,  ...,  on  the  interval  [0,T], 
where  h  =  -  t^  is  the  grid  spacing.  One  then  approximates  the 

~  nation  y(t)  of  Equation  2.26  at  the  time  grid  points  as  y1  «  y(t^). 

The  basic  approximating  equation  for  stiff  differential 
equations  [35  and  36]  is  the  Gear  algorithm 
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n+1 

y 


*0y 


n+1 


k 

I 

j=l 


n-j+1 


(2.28) 


where  the  constants  8^  and  G  ,  j=l,...,k,  called  Gear's  coefficients, 
are  determined  so  that  Equation  2.28  is  exact  for  any  polynomial 
solution  of  Equation  2.26  of  degree  up  to  k.  Gear's  coefficients 
[  35  and  36]  also  have  the  property  that  the  algorithm  tends  to  be 
stable,  even  for  stiff  differential  equations;  i.e.,  differential 
equations  with  widely  split  eigenvalues. 

A  multistep  formula  that  is  used  to  solve  Equation  2.27  is 
derived  from  Equation  2.28.  One  progresses  from  t^  to  t^+^  by 
solving  Equation  2.28,  together  with 


g(y 


n+1 

f 


.n+1 

y  > 


=  0 


(2.29) 


Using  the  Newton  formula  to  solve  Equations  2.28  and  2.29,  simul¬ 
taneously,  leads  to 


Ay 


(m) 


(m)  . (m) 

*  H 


(2.30) 


where  Ay^  =  -  y^  and  m  represents  the  iteration  number. 

The  time-step  counter  has  been  dropped  here  for  notational  simplicity. 
Substitution  of  Equation  2.30  into  Equation  2.28,  noting  that  the 
summation  term  of  Equation  2.28  remains  constant  at  each  iteration, 
yields  the  corrector  formulas 
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It  is  interesting  to  note  that  every  element  of  y  is  obtained  at  each 
time  step,  even  though  many  elements  of  y  may  not  appear  in  any  of 
the  equations • 

Since  Gear's  algorithm  uses  the  Nordsieck  vector  representation, 
the  predicted  Nordsieck  vector  at  the  present  time  step  is  obtained 
by  multiplying  the  Pascal  triangle  matrix  times  the  Nordsieck  vector 
at  the  previous  time  step  [35l-  Details  of  convergence  criteria, 
error  control,  and  procedures  for  this  technique  are  discussed  in 
references  20,  35,  and  36. 

(b)  Sparse  Matrix  Algebra 

The  system  of  nonlinear  algebraic  and  differential  equations 
defined  in  the  previous  sections  is  loosely  coupled.  For  two  reasons, 
however,  no  attempt  is  made  to  obtain  a  smaller  system  of  equations 
by  reducing  the  number  of  equations  through  elimination  of  general¬ 
ized  coordinates.  First,  the  Jacobian  matrix  formed  by  linearization 
of  the  coupled  equations  and  used  in  iterative  solution  by  Newton's 
method  is  sparse  and  can  be  very  efficiently  stored  and  decomposed. 
Secondly,  the  repetitive  nature  of  the  loosely  coupled  equations  re¬ 
sults  in  compact  and  efficient  computer  routines  for  evaluating  the 
nonlinear  equations  and  nonzero  matrix  entries.  Recently  developed 
sparse  matrix  algorithms  [37]  make  both  of  these  operations  prac¬ 
tical  and  desirable.  It  has  been  shown  [37]  that  is  is  usually 
more  efficient  to  solve  large  systems  of  sparse  equations ,  rather 
than  smaller  systems  with  a  greater  percentage  of  nonzero  entries. 


36 


Consideration  of  matrix  sparsity  is  important  to  the  speed  of 
computation  in  problems  of  dynamic  system  analysis  [35]-  when  less 
than  30/c  of  the  matrix  entries  are  nonzero,  it  is  inefficient  to 
store  the  matrix  as  a  two-dimensional  array.  Instead,  the  nonzero 
entries  are  stored  in  compacted  form.  The  method  most  commonly  used 
for  compacting  the  data  is  to  store  the  row  and  column  indices  of 
each  nonzero -valued  entry  in  the  matrix  in  two  vectors  i  and  j,  and 
its  value  in  a  third  vector  A.  This  is  called  "i-j"  ordering  and  is 
a  method  of  initially  storing  data  from  a  user-supplied  list  of 
physical  system  elements. 

Sparse  matrix  algorithms  are  most  efficient  when  the  nonzero 
matrix  entries  are  stored  in  an  organized  manner.  This  usually  im¬ 
plies  that  they  are  evaluated  row-by-row  or  column-by-column.  The 
previously  mentioned  repetitive  matrix  evaluation  scheme  is  defeated 
by  this  requirement,  since  it  usually  results  in  the  evaluation  of 
small  submatrices  located  at  various  positions  throughout  the  matrix. 
To  overcome  this  difficulty,  a  special  permutation  vector  is  generated 
from  the  row  and  column  vectors  describing  the  non zero -valued  posi¬ 
tions.  As  each  matrix  entry  is  generated,  it  is  directed  to  a 
specific  location  in  the  "A"  vector  by  a  permutation  index.  At 
completion  of  the  matrix  evaluation,  all  entries  are  stored  exactly 
as  if  they  had  been  evaluated  in  column  order.  A  sparse  matrix 
decomposition  algorithm  is  then  applied  to  the  column -ordered  matrix 
and  the  standard  LU  factorization  [35]  is  accomplished.  Full  pivot¬ 
ing  is  not  achieved,  but  the  algorithm  chooses,  among  the  acceptable 
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pivot  elements,  the  pivot  that  results  in  a  minimum  number  of  fills 
in  the  resulting  L  and  U  matrices.  This  is  important  for  efficient 
execution  of  the  forward  and  backward  substitution  phases,  since  an 
increased  number  of  fills  destroys  the  original  matrix  sparsity  and 
results  in  excessive  computer  time. 


■\K' 


CHAPTER  3 


KINEMATIC  ANALYSIS  AND  GENERALIZED  COORDINATE 
PARTITIONING  FOR  DIMENSION  REDUCTION  IN 
ANALYSIS  OF  CONSTRAINED  RIGID  BODY  SYSTEMS 

3.1  Introduction 

The  implicit  meth  1  of  simultaneously  solving  systems  of  differ¬ 
ential  and  algebraic  equations  of  motion  described  in  Chapter  2  is 
not  desirable  for  several  reasons.  Implicit  methods  employ  iterative 
techniques  requiring  a  Jacobian  matrix  of  the  combined  system  of 
Equations  2.8,  2.22,  and  2-23.  For  large  systems,  analytical  expres¬ 
sions  for  derivatives  of  Equation  2.23  are  required  if  program 
efficiency  is  to  be  achieved.  Generalized  forces,  possibly  discon¬ 
tinuous,  may  be  represented  by  digitized  data,  from  which  it  is 
difficult  or  impossible  to  provide  derivatives,  either  analytically 
or  numerically. 

In  the  implicit  method,  the  coordinates  q,  velocities  w,  and 
Lagrange  multipliers  \  of  Equations  2.22  and  2.23  are  treated  as 
state  variables,  which  are  determined  iteratively  by  the  integration 
algorithm.  The  algorithm  obtains  q  and  w  by  integration,  allowing 
it  to  maintain  error  control  and  accurate  prediction  of  these  vari¬ 
ables.  The  algorithm  uses  the  same  predictor  for  Lagrange  multipliers. 
Since  generalized  forces  may  be  highly  irregular  or  discontinuous, 
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irregular  and  discontinuous  joint  reaction  forces  are  reflected  in 
tile  Tiagrango  multipliers.  The  net  result  is  a  poor  prediction  of 
Lagrange  multipliers,  requiring  more  iterations  to  achieve  conver¬ 
gence  (or  quite  often  divergence)  in  the  corrector  step.  A  reduction 
in  time  stepsize  is  then  required  to  achieve  better  predicted  values. 

Stepsize  reduction  in  the  implicit  algorithm  is  undesirable 
since  it  requires  more  computer  time  and  often  leads  to  failure  of 
the  algorithm  to  achieve  a  solution,  even  when  extended  precision 
computer  arithmetic  is  used.  The  reason  for  corrector  divergence  is 

as  follows :  The  time  step  h  appears  in  the  denominator  of  the 

/  \ 

expression  multiplying  the  term  in  Equation  2.33,  and  as  h  gets 

smaller,  these  terms  get  larger  and  dominate  the  Jacobian  matrix. 

The  Jacobian  matrix  may  thus  become  badly  conditioned  or  singular. 

A  badly  conditioned  matrix  results  in  erroneous  Newton  differences 
associated  with  Lagrange  multipliers,  leading  to  divergence  of  the 
corrector. 

The  numerical  method  developed  in  the  following  chapters  offers 
significant  improvement  over  the  implicit  method  of  Chapter  2. 
Iterations  are  limited  to  the  solution  of  constraint  equations,  thus 
requiring  a  much  smaller  Jacobian  matrix.  Since  iteration  only 
involves  predicted  variables  with  error  control,  convergence  is  much 
faster.  In  fact,  the  iteration  step  is  frequently  bypassed,  because 
the  constraint  equations  are  satisfied  by  the  predicted  variables. 
Badly  conditioned  matrices,  as  described  earlier,  are  avoided  and 
single  precision  computer  arithmetic  can  be  used. 


Another  advantage  of  the  proposed  method  is  that  a  minimal  set 
of  differential  equations  is  identified  by  the  program  and  solved  by 
an  explicit  integration  method.  This  avoids  a  Jacobian  matrix 
involving  derivatives  of  generalized  forces. 

3-2  System  Constraints 

(a)  Classification  of  System  Coordinates 
The  ultimate  goal  of  this  research  is  to  develop  analysis 
techniques  that  use  system  constraints  to  improve  numerical  effi¬ 
ciency  and  accuracy,  while  at  the  same  time  reducing  the  preprocessing 
requirements  of  the  user.  To  this  end,  an  understanding  of  what 
constitutes  a  constraint  is  essential.  According  to  Webster  [38]  a 
constraint  is  confinement;  restriction;  force;  compulsion;  or 
coercion.  In  the  tradition  of  mechanics,  a  constraint  is  any  condi¬ 
tion  that  reduces  the  freedom  of  a  system.  Another  notion  is 
"something  that  limits  motion". 

In  the  development  of  Chapter  2,  a  system  was  assumed  to  be 
composed  of  n  rigid  bodies  in  a  plane,  each  with  three  degrees  of 
freedom.  Algebraic  constraint  equations  were  also  formulated  and 
reaction  forces  at  the  constraint  surfaces  were  introduced  into  the 
system  differential  equations  of  motion,  through  a  vector  of  Lagrange 
multipliers.  The  equations  of  motion  and  constraint  equations  were 
solved  iteratively  and  no  account  was  taken  of  the  number  of  system 
degrees  of  freedom.  Mathematically,  this  system  has  3n  degrees  of 
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freedom  and  the  algebraic  constraints/constraint  reaction  force 
system  represents  a  restriction  on  admissible  movements. 

The  formulation  of  Chapter  2  treats  every  coordinate  as  inde¬ 
pendent,  hence  the  system  has  3n  generalized  coordinates.  A 
formulation  reduces  the  number  of  degrees  of  freedom  when  the  alge¬ 
braic  constraints/constraint  reaction  force  balance  is  used  to 
eliminate  generalized  coordinates  from  the  system  equations  of  motion. 
Application  of  constraint  relations  to  mathematically  reduce  the 
number  of  degrees  of  freedom  introduces  a  new  problem  of  identifying 
independent  and  dependent  variables.  If  constraint  equations  mathe¬ 
matically  eliminate  r  degrees  of  freedom,  there  are  (3n-r)  independent 
generalized  coordinates  and  the  remaining  r  coordinates  become 
dependent  position  coordinates.  For  notational  purposes,  let  the 

symbol  u  designate  dependent  variables  and  let  v  designate  independent 

T  T  T 

generalized  coordinates.  With  this  notation,  q  =  [u  ,  v  ]  . 

Partitioning  of  q  into  u  and  v  must  be  done  in  a  manner  that 
does  not  degrade  program  efficiency  or  increase  numerical  errors.  If 
the  partitioning  were  done  on  a  random  basis,  there  are  3n!/rl (3n-r) I 
possible  combinations,  most  of  which  would  be  mathematically  unaccept¬ 
able.  Wells  [39]  stated,  when  discussing  the  selection  of 
independent  coordinates  for  simple  problems:  "it  is  a  well  known 
fact  that  certain  coordinates  may  be  more  suitable  than  others. 

Hence  the  quantities  chosen  in  any  particular  case  are  those  which 
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appear  to  be  the  most  advantageous  for  the  problem  in  hand.  The  final 
choice  depends  largely  on  insight  and  experience."  This  approach  is 
also  taken  in  another  well  developed  and  documented  planar  rigid  body 
modelling  program  [l3 ].  However,  in  this  case,  steps  are  taken  to 
circumvent  some  of  the  numerical  difficulties  associated  with  bad 
choices  of  independent  variables.  When  systems  become  large  and 
complex,  geometrical  insight  and  experience  quite  often  fail  and 
mathematical  techniques  must  be  employed. 

To  the  author's  knowledge,  the  first  organized  method  for 
selection  of  a  good,  and  possibly  best,  set  of  independent  generalized 
coordinates  for  large  scale  systems  was  developed  by  Sheth  [40,  l6] . 
Although  he  does  not  prove  that  the  selected  set  is  best,  he  does 
provide  numerous  arguments,  based  on  geometrical  considerations  and 
example  test  problems,  that  indicate  the  set  has  desirable  features 
of  a  properly  driven  kinematic  mechanism.  This  approach  has  sound 
mathematical  appeal  and  will  be  developed  as  the  tool  for  efficient 
dimension  reduction  of  constrained  systems. 

(b)  Iterative  Solution  of  System 
Constraint  Equations 

The  system  constraint  equations  developed  in  Chapter  2  are 
written  symbolically  in  the  form 


j(q,t)  =  0 


(5.1) 
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where  all  equations  may  not  be  independent,  but  they  are  consistent. 
That  is,  the  Jacobian  matrix  5^  with  m  rows  and  3n  columns  has  row 
rank  0  s  r  s  m  and  r  £  3n.  There  are  one  or  more  nonsingular  sub- 
matrices  of  s  of  rank  r.  Gauss -Jordan  reduction  of  the  matrix  $  , 

q  q 

T  T  T 

with  double  pivoting,  defines  a  partitioning  of  q  =  [u  ,  v  ]  such 

that  $u  is  a  submatrix  of  $  of  rank  r,  whose  columns  correspond  to 

elements  u  of  q,  and  is  a  submatrix  of  whose  columns  correspond 

to  elements  v  of  q.  Furthermore,  the  matrix  8  has  ideal  numerical 

u 

properties  associated  with  double  pivoting. 

The  subprogram  described  in  Appendix  A  determines  the  rank  of  ^ 
8  and  its  linearly  independent  rows  and  columns.  It  performs 
Gaussian  elimination  with  row  and  column  permutations  that  bring 
independent  rows  and  columns  to  the  top  and  left  of  the  matrix,  until 
the  remaining  submatrix  in  the  lower -right  corner  is  null  (all  of  its 
elements  are  less  than  some  specified  tolerance).  In  addition,  the 
matrix  is  decomposed  into  factors  of  the  form 


L 

0 

U 

UR 

y 

-  - 

0 

and 

0 

0 

which  are  superimposed  onto  the  original  matrix.  The  factor  products 


L-U 

L*UR 

_ 

LR  *U 

LR-UR 

represent,  in  permuted  form,  the  original  matrix  #  ;  specifically, 
L*U  =  #  and  L’UR  =  Jv>  The  columns  of  U  and  UR  determine  a 
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permutation  of  q  into  dependent  (u)  and  independent  (v)  variables, 

respectively.  The  rows  of  L  and  LR  determine  a  permutation  of  con- 

1  2 

straint  equations  §  =  0  into  independent  $  and  dependent  5  equations, 
respectively.  In  addition,  the  respective  lower  and  upper  triangular 
matrices  L  and  U  are  nonsingular. 

In  the  following  chapters  it  is  often  convenient,  for  clarity, 
to  express  matrix  equations  in  the  above  permuted  and  factored  forms, 
rather  than  in  terms  of  the  original  matrix  {  . 

Some  additional  comments  can  be  made  concerning  the  matrix  $ 

q 

and  its  corresponding  submatrices.  First  of  all,  may  be  null. 

That  is,  there  may  be  no  constraints  in  the  system,  in  which  case 
v  =  q  and  u  is  null.  This  special  case  requires  no  further  analysis 


and  need  not  be  considered  here. 


If  ?v  is  null  there  are  no  independent  generalized  coordinates 
and  u  =  q.  In  this  case,  in  view  of  Equation  3.1,  the  system  is 
either  a  structure  ?  =  ?(q)  or  its  geometry  is  completely  determined 
as  a  function  of  time,  $  =  $(q,t). 

If  roots  of  Equation  3.1  are  to  be  found,  an  iterative  procedure 
such  as  Newton’s  method  is  required.  With  an  initial  estimate  q(0) 
of  the  system  configuration  at  time  tQ,  the  matrix  ?  (q(0),tQ)  is 
evaluated.  The  iterative  equation  is  then 


«q[q(j),t0]  [q(j+l)  -  q  (<j )  ]  =  -  5[q(j),t0],  j  =  0,1,...  (3.2) 
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where  j  is  an  iteration  counter.  The  solution  of  Equation  3.1 
requires  the  determination  of  q(j+l)  in  Equation  3.2.  If  Equation 
3.1  is  consistent,  there  is  at  least  one  solution.  If  the  rank  is 
less  than  3n  there  are  an  infinite  number  of  solutions . 

Assume  that  the  matrix  has  been  decomposed  by  the  subprogram 
described  in  Appendix  A.  Equation  3*2  may  then  be  written  in  per¬ 
muted  form,  in  terms  of  the  previously  mentioned  matrix  factors  and 
partitioning,  as 

L  •  U  •  [u(j  +  1)  -  u(j)  ]  +  L  •  UR  •  [v(j  +  l)  -  v(j) ] 

=  *  (3.3) 

LR  ■  U  •  [u(j  +  1)  -  u( j) J  +  LR  •  UR  •  [v(j  +  l)  -  v(j) ] 

=  -  ?2[q(l),t0]  (3.4) 

Noting  that  v  are  independent  generalized  coordinates,  one  may  set 
v( j + l)  =  v(j)  and  Equation  3.3  reduces  to 

L  ■  U  •  [uU  +  1)  -  u(j)J  =  -  i1  (3.5) 

Since  matrices  L  and  U  are  nonsingular  and  triangular,  simple 
forward  elimination  and  back  substitution  determines  u(j  +  l). 
Equation  3.4  is  not  required  in  the  intermediate  steps  of 
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determining  roots  of  Equation  3.1,  but  may  be  used  to  test  for  con- 
sistency  of  constraints  at  the  final  step.  That  is,  ||$  \\  must  be 
less  than  a  specified  tolerance,  if  redundant  constraints  are  to  be 
satisfied. 

Kinematic  analysis  of  constrained  systems  requires  that  bodies 

move  consistent  with  constraints.  If  the  system  is  moved  from  posi- 

1  2 
tion  q  at  tq  to  a  nearby  position  q  at  t^,  consistent  with 

X  2  2 

constraints,  then  f(q  ,tq)  =  $(q  ,t^)  =  0.  If  §(q  ,tg)  is  expanded 
in  a  Taylor  series  about  q1  and  t  one  has 

$(q2,t2)  =  fCq1,^)  +  $qAq  +  $tAt  +  |  (*qAq)q  Aq 

+  ...  =  0  (3.6) 

2  1 

where  Aq  =  q  -  q  and  At  =  -  tq.  If  Aq  and  At  are  small,  second 

and  higher  order  terms  of  Equation  3.6  can  be  neglected,  yielding  a 

2 

linear  predictor  of  q  ,  of  the  form 

$(q2  -  q1)  =  -  «+At  (3.7) 

q  ^ 

Comparing  Equations  3. 7  and  3.2  and  noting  the  permuted  form  of 
Equations  3.3  and  3.4,  one  can  express  Equation  3.7  as 

L  •  U  •  (u2  -  u1)  +  L  •  UR  •  (v2  -  v1)  =  -  }*  At  (3.8) 


hi 

LR  •  U  •  (u2  -  u1)  +  LB  •  UR  •  (v2  -  v1)  =  -  $2  At  (5-9) 

2 

Since  L  and  U  are  nonsingular  and  with  v  and  t£  specified.  Equation 
2 

3.8  determines  u  .  It  is  most  efficient  to  solve  Equation  3.8  in  two 
steps,  since  the  matrix  product  L  •  UR  is  not  calculated  by  the  sub¬ 
program.  First,  solve 

L  •  U  •  u*  =  -  § l  At  (3.10) 

Tj 

by  forward  elimination  and  back  substitution  and  then  use  the  pre¬ 
calculated  matrix  H  to  obtain 

u2  -  u1  =  H-(v2  -  v1)  +  u*  (3.11) 

2 

Equation  3.3  is  then  used  to  iteratively  correct  u  ,  to  satisfy 
Equation  3.1. 

The  matrix  product  L  *  UR  of  Equation  3.8  is  never  evaluated, 
since  the  matrix 

H  =  -  U-1  •  L_1  *  L  •  UR  =  -  U_1  •  UR 

is  required  in  Equation  3-1.  In  fact  all  subsequent  equations  in¬ 
volving  the  matrix  L  •  UR  ultimately  require  the  matrix  H  for 
solution,  so  for  convenience,  H  is  calculated  and  stored  in  the  space 
originally  occupied  by  UR. 
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To  avoid  excessive  corrector  iterations,  it  is  desirable  that 
2 

Equation  3.1  predict  u  with  small  error.  Since  Aq  was  assumed  small 
and  Av  is  fixed,  Au  should  be  small.  Gaussian  elimination,  using 
full  row  and  column  pivoting  to  partition  q  into  u  and  v,  accomplishes 
this.  Equation  3.1  may  be  written  simply  as 

* 

Au  =  H  Av  +  u 


As  noted  in  Appendix  A,  the  largest  matrix  element  is  permuted  to  the 
diagonal  position  at  each  intermediate  step,  when  performing  Gaussian 
elimination  with  full  row  and  column  pivoting.  This  process  tends  to 
maximize  the  magnitude  of  diagonal  elements  of  U  and  to  minimize  the 
magnitude  of  elements  of  UR.  Noting  that  H  =  -  U  1  *  UR,  one  would 
expect  that  the  magnitude  of  elements  in  H  will  be  substantially  less 
than  1  and  that  the  magnitude  of  elements  in  Au  will  be  substantially 
less  than  the  magnitude  of  elements  in  Av. 

This  argument  could  be  strengthened  if  a  reliable  estimate  for 
the  condition  number  of  U,  cond(U) ,  were  available.  In  general, 


Gaussian  elimination  with  full  row  and  column  pivoting  yields  well 
conditioned  matrices  [4lJ .  However,  one  can  find  examples  in  which 
it  does  not.  If  cond(u)  and  the  norm  of  U  are  known,  the  norm  of 


cond(u)/| |u| | 


(3.1?) 


* 


U9 


The  norm  of  U  is,  by  construction,  large  since  the  largest  elements  of 
5  were  permuted  into  it.  If  cond(u)  is  small  (on  the  order  of  l) , 
then  by  Equation  3-12,  | |u  1|  j  «  l/| |uj | . 

Since  the  elements  of  U  dominate  those  in  UR,  it  is  reasonable 
to  expect  that  | |u| |  »  | | UR [ |  and  | | U-1 | |  •  | | UR | |  «  1,  by  the 
assumption  of  a  well  conditioned  U  matrix.  Finally, 

| 1 H | j  =  |1-  U  1  •  UR | |  <  | |U  X| |  •  | | UR | |  «  1,  supports  the  contention 
that  | | Aq | |  is  reduced  by  proper  selection  of  u  and  v. 

3.3  Kinematic  Velocity  and 
Acceleration  Analysis 

(a)  Velocity  Analysis 

For  kinematic  analysis,  if  one  specifies  independent  general¬ 
ized  coordinates  as  functions  of  time,  independent  velocities  and 
accelerations  can  also  be  expressed  as  functions  of  time. 

The  remaining  dependent  velocities  and  accelerations  can  then 
be  determined  by  differentiating  the  constraint  equations  with  respect 
to  time.  Differentiating  Equation  3.1  with  respect  to  time  yields 

Sq  4  +  »t  =  0  (5.13) 

Comparing  Equations  3.13  and  3.2,  equations  similar  to  Equations  3-3 
and  3-4  can  be  expressed  as 


L*U*u  +  L*UR*v=-?^ 


(3.14) 


I 
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LR*U*u+LR-UR*V=-f^  (3.15) 

The  solution  of  Equation  3.1*+  determines  u  and  Equation  3.15  is  auto¬ 
matically  satisfied  by  the  consistency  requirement  of  Equation  3 .4. 

Similar  to  the  method  of  solving  Equation  3.8,  first  solve 

L  •  U  •  u*  =  -  (3.16) 

by  forward  elimination  and  back  substitution  and  then  use  the  pre¬ 
calculated  matrix  H  to  obtain 

u  =  H  v  +  u  (3.17) 

Inspection  of  Equation  3.17  reveals  that  the  matrix  H  relates 
dependent  velocities  to  independent  velocities  and  this  matrix  is 
generally  referred  to  as  an  influence  coefficient  matrix.  From  a 
geometrical  point  of  view,  a  properly  driven  mechanical  system  is  one 
in  which  the  ratios  of  output  to  input  velocities  are  minimized 
[40,  l6]  .  In  addition,  the  system  is  driven  with  the  greatest 
mechanical  advantage. 

Recall  that  the  method  of  decomposing  determined  H,  such 
that  Au  was  reduced  in  Equation  3.11,  hence  u  is  also  reduced  in 
Equation  3.17.  A  further  advantage  of  minimizing  velocity  ratios  is 
that  errors  in  independent  velocities  are  not  amplified  when  calcu¬ 
lating  dependent  velocities.  This  will  be  important  when  solving  the 
system  differential  equations  of  motion  in  Chapter  4. 
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(b)  Acceleration  Analysis 

Having  determined  all  system  position  and  velocity  variables, 
one  may  proceed  to  determine  dependent  accelerations.  Equation  5.13 
may  be  differentiated  with  respect  to  time  to  obtain 

H+  (8  d)  q  +  2*  . q  +  I..  =  0  (3-18) 

q  q  q  qt)  bt 

For  convenience,  denote  the  three  known  terms  in  Equation  3.18 
as 

C(q,q,t)  =  (8qq)q  q  +  25qtq  +  5tt  (3.19) 

and  partition  Equation  3.18  as 

L-U*U+L*UR*v=-C1  (3.20) 

LR*U-u+LR'UR-v=-C2  (3.2l) 

The  independent  generalized  accelerations  v  are  presumed  known,  so 
Equation  3.20  determines  u  and  Equation  3.21  is  automatically  satis¬ 
fied  by  the  consistency  requirement  of  Equation  3.4. 

Similar  to  Equation  3.8,  Equation  3.20  is  solved  in  two  steps 

*  1 
L  *  U  '  u  =  -  C 


(3.22) 
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by  forward  elimination  and  back  substitution.  Also, 


* 

u  =  H  v  +  u 


(3.23) 


where  the  precalculated  influence  coefficient  matrix  H  is  used.  Note 
that  errors  in  independent  generalized  accelerations  v  are  reduced, 
in  calculating  the  dependent  accelerations  h,  because  of  the  small 
norm  of  H. 


3.4  Kinetostatics  of 
Constrained  Systems 

The  analysis  method  of  Section  3.2  allows  one  to  identify  an 
ideal  set  of  independent  generalized  coordinates,  from  which  to  drive 
an  arbitrary,  constrained  multiple -degree -of -freedom  system.  Kineto- 
static  analysis  may  then  be  applied,  using  the  method  of  Section  3-3, 
to  determine  all  forces  in  the  system,  including  inertial  forces, 
generalized  joint  reaction  forces,  and  the  generalized  driving  forces 
required  to  achieve  the  specified  system  dynamics.  One  could  then 
envision  the  same  system  reacting  to  this  set  of  generalized  driving 
forces.  In  this  case,  a  system  of  differential  equations  of  motion, 
solved  by  the  methods  of  Chapter  4  will  determine  the  same  dynamic 
response.  The  first  system  is  displacement  driven,  requiring  speci¬ 
fied  generalized  coordinates,  and  the  second  system  is  force  driven, 
requiring  specified  generalized  forces.  The  two  examples  demonstrate 
a  duality  between  specified  generalized  coordinates  and  specified 


generalized  forces. 
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Using  the  permutations  of  the  previous  sections,  equations  of 
motion  for  constrained  systems  derived  in  Chapter  2  may  be  partitioned 
into 


1  T 
X  +  U 


,  2  ,,uu..  ,  „„uv  .. 
X  =  M  u  +  M  V 


-  (AU  +  QU)  (3.2k) 

URT  •  LT  •  X1  +  URT  •  LRT  •  X^  =  +  MWV- 

-  (AV  +  QV)  (3.25) 

1  2 

Lagrange  multipliers  X  and  X  are  associated  with  independent  con¬ 
straint  equations  =  0  and  redundant  or  dependent  constraint 

2  2 
equations  5  ,  respectively.  The  multipliers  X  are  not  determined  by 

Equations  3.2k  and  3.25,  but  must  be  obtained  from  other  methods, 

such  as  local  body  deformation  models  that  specify  relations  between 

X^  and  X2. 

Let 

FU  =  U1  •  LT  ’  X1  +  UT  •  LRT  •  X2  (3.26) 


(3.27) 


Then,  by  combining  Equations  3.26  and  3.27  with  H  =  -  U  1  •  UR,  it  is 
easy  to  show  that 


v  T  u 

F  =  -  H  F 


(3.28) 


1  2 

holds  for  any  value  of  \  and  \  .  Assuming  that  redundant  constraints 

are  consistent,  they  can  be  eliminated  without  effecting  the  system 

2 

dynamics,  which  is  then  equivalent  to  setting  \  to  zero  in  Equations 
3 • 2k  and  3  •  25  • 

A  displacement  driven  system  with  specified  generalized  coor¬ 
dinates  has  corresponding  unknown  generalized  forces.  In  this  case, 

v  v* 

let  the  generalized  force  vector  Q  be  composed  of  known  Q  and 
unknown  RV  forces 


Q 


V 


+  R 


(3.29) 


Generalized  constraint  reaction  forces  become 


FU  =  UT  •  LT  •  A1  =  MUU\1  +  MUVv  -  (AU  +  QU) 


(3.30) 


and  the  unknown  generalized  driving  forces  become 
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It  is  not  necessary  to  evaluate  \  to  determine  the  generalized  con¬ 
straint  reaction  forces . 

Equations  3.3 0  and  3.31  may  be  made  more  tractable  by  substi¬ 
tuting  Equation  3.23  to  eliminate  dependent  accelerations,  thus 

FU  =  (MUUH  +  MUV)V  -  (AU  +  QU  -  MUUu*)  (3-32) 

RV  =  (MVUH  +  -  (AV  +  QV*  -  M^*)  +  HTFU  (3.33) 

If  forces  FU  are  not  desired,  one  may  combine  Equations  3.32  and  3.33 
to  eliminate  FU,  obtaining 

sv  =  [Mw  +  M™H  +  hVV  *  MuuH)]v  -  (AV  +  QV*  -  m'™!!*) 

-  HT(Au  *  «u  -  MUV)  (3.5U) 

Equation  3.3*+  thus  relates  generalized  driving  forces  RV  to  specified 
generalized  coordinates  v,  where  velocities  and  accelerations  are 
assumed  known. 

The  kinetostatics  problem  with  full  consideration  of  all 
forces  in  the  system  may  now  be  solved.  The  basic  steps  are,  pre¬ 
suming  equations  and  initial  coordinate  estimates  have  been  given: 

Step  1.  Identify  dependent  and  independent  variables. 
Decompose  the  matrix  8  ,  using  the  method  of  Appendix  A,  to  determine 
a  partitioning  of  q  into  u  and  v. 
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Step  2.  Solve  the  constraint  equations.  Specify  v  and  iter¬ 
atively  solve  Equation  3.1,  using  Equation  3. 5*  Check  consistency  of 
2 

equations  $  =0. 

Step  3.  Calculate  system  velocities.  Specify  v  and  determine 
u  from  Equations  3-l6  and  3.17- 

Step  4.  Calculate  system  accelerations.  Specify  V,  evaluate 
C  ,  and  determine  u  from  Equations  3.22  and  3.23.  Dependent  acceler¬ 
ations  U  are  not  required  if  Equations  3.32,  3.33,  or  3*3^+  are  used 
in  Step  5. 

Step  5.  Determine  generalized  constraint  reaction  forces  FU 

v  u 

and  system  driving  forces  R  .  Evaluate  A  and  Q,  and  determine  F  from 

Equation  3.32  and  RV  from  Equation  3.33  or  Equation  3.3*+.  Individual 

constraint  reaction  forces  are  obtained  by  first  solving 


T  T  1 
U  L  \  = 


(obtained  from  Equation  3.26  with  \ 


2 


0)  for  k  .  Note  that 


F 


u 


=  5 


V 

u 


■v 

V 


and 
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T  T 

Each  column  entry  in  a  given  row  of  ($v)  times  the  corresponding 

row  entry  of  A.1  determines  the  generalized  constraint  reaction  force 

contribution  to  the  generalized  coordinate  associated  with  the  given 
T  T 

row  of  (5  ).  Hence,  one  may  calculate  reaction  forces  at  any  con¬ 
straint  in  the  system. 

5.5  Sparse  Matrix  Considerations 
As  noted  in  Chapter  2  the  nodal  method  of  constraint  formula¬ 
tion  yields  loosely  coupled  equations  with  a  corresponding  sparse 
Jacobian  matrix.  If  the  system  has  a  large  number  of  constraints, 
sparse  matrix  manipulation  algorithms  will  be  efficient.  While  the 
equations  derived  in  previous  sections  still  apply,  they  may  be 
expressed  in  different  forms  to  avoid  unnecessary  calculation  of 
sparse  matrix  products. 

Sparse  matrix  algorithms  that  factor  matrices  are  not  suitable 
for  determining  the  partitioning  of  coordinates  q  into  u  and  v.  They 
employ  a  partial  pivoting  strategy  and  to  maintain  efficiency,  usually 
do  not  select  the  largest  pivotal  elements.  The  matrix  so  identi¬ 
fied  will  have  a  smaller  norm  than  the  corresponding  matrix  identified 
by  full  pivoting.  The  corresponding  matrix  H  =  -  will  then 

have  a  larger  norm.  Therefore,  Gaussian  elimination  with  full  pivot¬ 
ing  is  periodically  applied  to  the  full  matrix  representation  of 
solely  for  the  purpose  of  determining  the  ordering  of  q  into  u  and  v. 
Once  established,  a  matrix  I  is  appended  to  f  to  form  a  modified 
matrix  $  of  rank  5n  such  that  when  $  is  permuted,  corresponding  to 

1  q 
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the  permutation  of  into  and  it  is  of  the  form 


where  I  is  an  identity  matrix.  The  matrix  I  is  readily  constructed 
from  the  permutation  vectors  generated  by  Lhe  algorithm  described  in 
Appendix  A. 

An  algorithm  employing  sparse  matrix  techniques  for  a  dis¬ 
placement  driven  system  with  full  consideration  of  inertial  and 
externally  applied  loads  is  given  in  the  following  steps: 

Step  1.  Periodically  identify  dependent  and  independent 
variables  by  decomposing  the  matrix  $  ,  using  the  method  of  Appendix 
A.  Construct  the  matrix  by  appending  the  matrix  I  to  {^. 

Step  2.  If  a  first  order  prediction  of  dependent  variables 
1  2 

(u  )  to  an  adjacent  position  (u  )  is  desired,  append  to  Equation  5.7, 

5  Aq  =  -  $,  At 
q  t 

the  equation 

.  2  1 

Av  =  v  -  v 

which  in  matrix  form  becomes 


(5.55) 


The  matrix  °f  rank  r  =  3n  is  then  factored  by  sparse  matrix 
algorithms  and  Aq  is  determined  by  Equation  3.35* 

Step  3.  The  predicted  dependent  variables  u  are  now  corrected 
iteratively  until  the  constraint  equations  are  satisfied  by  appending 
to  Eqiiation  3.2, 

$  Aq  =  -  § 

q 

the  requirement  that  independent  variables  remain  fixed 

Av  =  0 

which  in  matrix  form  becomes 

V<!-['o]  <5'36> 

Dependent  variables  are  iteratively  corrected  until  the  independent 
constrinats  are  satisfied  and  any  dependent  constraint  equations  are 
then  checked  for  consistency. 

Step  k.  Analysis  of  displacement  driven  systems  requires  that 

•  s 

independent  velocities  be  specified.  Denoting  these  as  v  ,  write 
Equation  3.13  as 


6o 


and  append  the  equations 


•  -s 
v  =  v 


which  in  matrix  form  becomes 


(3.57) 


Equation  3*37  then  determines  q. 

Step  5.  Combine  Equation  3- 18  and  v  =  vS  to  form  the  matrix 
equation 


(3.38) 


where  C  (q,  q,  t)  is  given  by  Equation  3.19*  Solve  for  accelerations 

q* 

Step  6.  Using  Equation  2.23  and  noting  from  Equation  3.29 

y 

that  unknown  generalized  forces  R  associated  with  independent 
generalized  coordinates  must  be  determined,  one  has 


5^  \  +  Is  Pv  =  Mq  -  (A  +  Q) 


which  in  matrix  form  becomes 


=  Mq  -  (A  +  Q) 


r  ^9) 


6l 


As  noted  in  Section  3.^,  Lagrange  multipliers  corresponding  to 
redundant  constraints,  are  not  determined  by  Equation  3.39>  However, 
if  redundant  constraints  are  consistent,  their  corresponding  multi¬ 
pliers  can  be  set  to  zero.  The  remaining  multipliers  and  unknown 
forces  RV  are  then  determined  by  Equation  3-39  and  all  forces  in  the 
system  are  determined. 

Equations  3-35  to  3.39  employ  the  same  sparse  matrix  5  or 

— T 

5  ,  thus  only  one  matrix  factorization  is  required  for  each  execution 
of  steps  2  to  6.  Matrix  and  vector  products  are  also  indicated  for 
evaluation  of  the  right  hand  side  of  Equations  3.38  and  3.39.  How¬ 
ever,  these  matrices  are  very  sparse  and  analytical  expressions  are 
written  directly  for  the  indicated  products  requiring  no  additional 
matrix  manipulation. 

To  maintain  program  efficiency  the  full  matrix  decomposition 
of  Step  1  should  not  be  performed  too  frequently,  yet  often  enough 
to  avoid  numerical  problems  associated  with  nearly  singular  5^ 
matrices.  Available  methods  for  monitoring  the  condition  number  of 
5^  require  more  time  than  is  required  for  its  original  determination. 
Monitoring  velocity  ratios,  however,  may  be  an  inexpensive,  though 
not  entirely  reliable,  technique  for  determining  the  need  for  a  new 
5^  matrix.  Equation  3.13  may  be  partitioned  and  written  as 

•  -1  • 

u  =  -  s  fa  v+5i 
L*v  *t. 
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To  approximate  the  amplifying  ability  of  one  could  evaluate  the 
number 


following  the  initial  determination  of  $  ,  (i  =  0)  and  at  each 
successive  application  of  the  sparse  matrix  code  (i  >  0) .  When  the 


estimate  of  increase  in  matrix  norm  K^/Kq  exceeds  a  specified  level. 


a  need  for  a  new  partitioning  of  5 


q 


is  then  indicated. 
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CHAPTER  4 

CONSTRAINED  EQUILIBRIUM  AND  DYNAMIC  ANALYSIS 

4.1  Introduction 

In  Chapter  3,  differential  equations  of  motion  axe  developed 
for  a  general  constrained  planar  rigid  body  system.  It  was  demon¬ 
strated  that  the  complete  system  state  is  determined  when  all 
independent  generalized  coordinates,  velocities,  and  acclerations 
are  specific  functions  of  time.  These  equations,  however,  must  be 
numerically  integrated  when  independent  variables  are  allowed  to  move 
under  the  influence  of  applied  force.  Therefore,  this  chapter  is 
devoted  to  development  of  an  efficient  numerical  integration 
algorithm  utilizing  kinematic  properties  of  constraints  to  obtain 
a  minimal  system  of  differential  equations  of  motion. 

Efficiency  is  achieved  by  minimizing  the  number  of  differential 
equations  to  be  solved,  which  improves  numerical  integration  perfor¬ 
mance.  Paul  [42]  has  indicated  that  the  minimum  number  of  first 
order  differential  equations  formulated  by  any  of  the  existing 
general  purpose  programs  is  2  dof  (degree  of  freedom) .  Only  2  dof 
first  order  equations  are  required  by  the  method  described  in  this 


chapter . 
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4.2  Equilibrium  Analysis 

Transient,  dynamic  analysis  of  complex  constrained  mechanical 
systems  is  often  initiated  from  a  position  of  partial  or  complete 
equilibrium.  Static  equilibrium  analysis  is  often  required  prior  to 
initiation  of  transient  analysis  of  vehicular  systems,  because  the 
nonlinear  characteristics  of  suspension  elements  make  it  difficult 
to  estimate  equilibrium  configurations.  Partial  equilibrium  analysis 
is  required  when  one  or  more  of  the  bodies  in  a  system  have  no 
equilibrium  position  or  when  internal  forcing  elements  such  as  springs, 
are  to  be  given  initial  lengths  differing  from  their  equilibrium 
lengths.  In  this  case,  temporary  constraints  are  placed  on  or 
between  bodies  prior  to  equilibrium  analysis  and  released  prior  to 
transient  analysis. 

Various  methods  may  be  employed  to  solve  equilibrium  equations . 
Since  they  are  highly  nonlinear,  iterative  techniques  are  required. 

One  iterative  technique  that  is  often  employed  is  Newton's  method 
for  finding  roots  of  nonlinear  algebraic  equations.  Newton's  method 
may  converge  to  unstable  equilibrium  configurations,  or  diverge,  if 
poor  initial  estimates  of  configuration  are  given.  However,  it  is 
easy  to  implement  and  has  been  used  to  find  equilibrium  of  systems 
with  many  degrees  of  freedom. 

The  implicit  numerical  integration  method  for  transient 
analysis  of  Chapter  2  requires  the  Jacobian  matrix  of  the  equations 
of  motion  for  solution  of  Equation  2.25.  Noting  that  velocities  and 
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accelerations  axe  zero  for  static  equilibrium.  Equation  2.23  expresses 
equilibrium  equations  as 


tQ)  X  +  Q(q,  tQ)  =  0 

> 

§(q,  tQ)  =  0 

/ 

and  the  corresponding  matrix  iterative  equation  is 


T  T 

(?  x  +  q)  $ 

q  q  q 

Aq(  j) 

T 

5q  X  +  Q, 

1  ** 

n 

o 

L 

AX(j) 

— 1 

5 

(4.1) 


(4.2) 


This  method,  requiring  the  simultaneous  determination  of  q 
and  X,  is  undesirable  because  it  involves  the  repetitive  solution  of 
a  large  system  of  equations  and  a  reasonable  estimate  of  the  vector 
X  must  be  given,  since  it  appears  in  the  matrix.  Poor  estimates  of 
X  may  lead  to  a  badly  conditioned  or  singular  matrix  in  Equation  4.2 
and  divergence  of  Newton's  method.  For  this  reason,  iterative 
techniques  that  do  not  involve  Lagrange  multipliers  are  desirable. 

In  Chapter  3  a  method  for  identifying  the  independent  gener¬ 


alized  coordinate  vector  v  was  developed  and  a  reduced  system  of 
differential  equations  of  motion,  Equation  3*34,  was  derived.  Noting 
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that  velocities  and  accelerations  are  zero  for  static  equilibrium, 
Equation  3-34  yields 

RV  =  -  H(q,  t0)T  Q(q,  tQ)U  -  Q(q,  t^^*  (4.3) 

where 

* (q>  tQ)  =  0  (4.4) 

v 

The  forces  R  must  be  applied  to  the  independent  or  specified 

g 

coordinates  v  in  a  displacement  driven  system.  In  an  equilibrium 

v 

situation  the  variables  v  are  free  and  the  forces  R  must  be  zero. 

A  condition  for  static  equilibrium  is  then 

RV  =  0  (4.5) 

Equations  4.4  and  4.5  may  be  solved  iteratively,  by  noting  that 
RV  =  RV(u,  v) .  Therefore,  with  du  =  Hdv,  the  iterative  equation 

f"  ^  +  H]Av  =  -  RV  (4.6) 

determines  the  solution  of  Equation  4.5.  The  matrix  in  Equation  4.6 
is  evaluated  by  numerical  differencing,  because  l)  the  matrix  H  and 


generalized  forces  appear  in  liquation  4.3,  2)  analytical  expressions 
for  their  partial  derivatives  are  difficult  to  obtain,  and  3)  the 
derivatives  are  not  required  for  subsequent  transient  analysis. 

4.3  Initial  Conditions 

Obtaining  suitable  initial  conditions  for  transient  analysis 
poses  problems  when  constrained  mechanical  systems  have  many  degrees 
of  freedom.  Often  one  desires  to  specify  initial  conditions  on 
various  coordinates  and  their  first  derivatives  or  initial  conditions 
on  first  derivatives  of  other  coordinates.  The  independent  variable 
set  identified  by  factoring  the  constraint  Jacobian  matrix  may  not 
agree  with  the  above  selected  variables  so  some  initial  conditions 
would  be  changed  when  satisfying  constraints.  The  method  described  in 
Section  3.6  allows  one  to  hold  selected  coordinates  and  velocities  at 
their  specified  values  by  appending  temporary  constraints.  The  only 
requirement  is  that  these  constraints  be  consistent  with  other  con¬ 
straints.  This  method  is  employed  below.  The  following  seven  steps 
outline  an  algorithm  for  achieving  suitable  initial  conditions,  from 
a  starting  configuration  that  is  consistent  with  constraints  or  from 
a  partial  or  ^ull  equilibrium  configuration: 

Step  1.  Specify  an  initial  estimate  q°  of  the  constrained 
system  configuration.  This  estimate  must  be  near  the  expected  final 
conf iguration  to  insure  convergence.  If  static  equilibrium  is  to  be 
obtained,  the  estimate  must  be  near  the  desired  equilibrium 


configuration  to  avoid  convergence  to  alternate  stable  or  unstable 
configurations,  or  divergence.  Set  q  =  q°. 

Step  2.  If  some  elements  of  q  are  to  remain  fixed  during 
initial  assembly  of  the  constrained  system  or  during  static  equilib¬ 
rium,  for  example  in  partial  equilibrium,  append  to  Equation  4.4 
additional  algebraic  constraints  of  the  form  q^  -  q^  =  0,  where  j 
identifies  the  variables  to  be  held  fixed.  Other  types  of  algebraic 
relations  between  coordinates  can  also  be  formulated  as  needed.  These 
constraints  will  be  removed  prior  to  step  6. 

Step  3.  Evaluate  the  modified  constraint  matrix  5  (including 

-  q 

constraints  of  step  2)  and  factor,  as  described  in  Appendix  A,  to 
determine  a  submatrix  of  maximal  rank  and  to  determine  a  partition¬ 
ing  of  q  into  u  and  v  (dependent  and  independent  coordinates, 
respectively) .  If  does  not  have  full  row  rank,  then  too  many 
constraints  or  inconsistent  constraints  are  specified.  Therefore 
some  equations  are  dependent  and  they  can  be  temporarily  ignored  in 
the  iterative  procedure  of  Step  4. 

Step  4.  Keeping  the  independent  variables  v  fixed,  iterate 
with  Equation  3.5  to  solve  the  independent  constraint  equations  of 
step  2,  until  1)?||  <  e,  where  e  is  a  specified  closure  tolerance 
level.  If  there  are  dependent  constraint  equations,  check  them  for 
violation.  When  dependent  constraints  are  violated,  inconsistent 
constraint  conditions  have  been  specified  and  a  meaningless  solution 
may  result.  If  static  equilibrium  is  not  desired,  go  to  step  6. 
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Step  5 .  Determine  a  new  estimate  for  v  from  Equation  4.6  that 
reduces  ||RV(|.  If  \ |Rvj  j  >  6,  a  specified  force  unbalance  tolerance, 
return  to  step  3.  If  the  iteration  in  steps  3  to  5  diverges,  either 
inconsistent  constraints  have  been  imposed  or  initial  estimates  of 

system  configuration  are  in  error. 

.  *0 
Step  6.  Specify  initial  estimates  for  velocities  q  and  set 

q  =  q^.  The  velocities  q  are  adjusted  to  satisfy  Equation  3.13, 

0 

evaluated  at  q  from  the  above  steps.  If  some  elements  of  q  are  to 

remain  fixed  at  q^,  append  to  Equation  3.13  additional  equations  of 
.  .  0 

the  form  q  -  q  =  0,  where  k  identifies  the  variables  to  be  held 

K  K 

fixed.  Other  types  of  relations  between  velocities  can  also  be  form¬ 
ulated  as  needed.  These  constraints  will  be  removed  prior  to  step  1 
of  transient  analysis. 

Step  7.  Factor  the  modified  matrix  of  Equation  3.13  as 
described  in  Appe  ,dix  A,  to  determine  a  partitioning  of  q  into  u  and 
v  (not  necessarily  the  same  partitioning  as  in  step  3).  Calculate  u 
using  the  above  partitioning  in  Equations  3.l6  and  3.17*  If  the 
above  matrix  lacks  full  row  rank,  some  equations  are  dependent  and 
must  be  checked  for  consistency.  If  dependent  equations  are  not 
satisfied,  inconsistent  constraints  on  velocity  are  specified  and  a 
meaningless  solution  may  result.  Otherwise,  a  consistent  set  of 
initial  conditions  is  now  available  for  transient  analysis. 

4.4  Numerical  Integration  of  the 
Equations  of  Motion 

The  reduced  system  of  differential  equations  of  motion, 
derived  in  Chapter  3,  were  expressed  in  the  form 
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RV  =[MW  +  M™H  +  HT(MUV  +  MUU  H)]v  -  (AV  +  -  M™  U*) 

-  HT(AU  +  QU  -  MUU  li*)  (U.7) 

Noting  that  the  independent  variables  v  are  free,  the  driving  forces 
RV  are  zero  in  Equation  U.7.  Therefore  Equation  U.7  can  be  written 
as 


[MW  +  HTMUV  +  (MW  +  HTMUU)H]V  =  AV  +  Q.V  +  HT(AU  +  QU) 

-  (M™  +  HTMUU)u*'  (U.8) 

In  order  to  integrate  Equation  U.8,  it  is  helpful  to  write  it 
in  first  order  form.  To  do  this,  define  the  vector  s  of  independe  b 
velocities  as 


v  =  s 


(4.9) 


and  write  the  total  vector  q  of  generalized  velocities  in  terms  of 
s,  using  Equations  3.16  and  3.17,  as 


w(s)  =  q(s)  = 


,* 

u 

Hs  +  u 

V 

s 

(4.10) 


.*  1 

where  L  •  U  •  u  «*  -  With  this  notation  Equation  3.19  becomes  a 
function  of  w, 
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C  =  ($w)w+2$  w  +  5+. 

q  q  qt  titi 

x- 

Equation  3.22  is  used  to  evaluate  u  ,  and  Equation  4.8  becomes 

[MW  +  HTMUV  +  (M™  +  HTMUU)H]s  =  AV  +  QV  +  HT(AU  +  QU) 

-  (M™  +  HTMUU)U*  (4.11) 

Decomposing  the  above  mass  matrix  as  described  in  Appendix  A,  one 
may  solve  Equation  4.11  for  s.  Equations  4.9  and  4.11  now  form  a 
system  of  6n  -  2r  =  2  dof  first  order  ordinary  differential  equations 
where  n  is  the  number  of  bodies  in  the  system  and  r  is  the  number  of 
dependent  variables  in  u. 

The  selection  of  numerical  integration  algorithms  for  dynamic 
analysis  should  be  based  on  the  type  of  differential  equations  to  be 
solved.  Since  most  mechanical  systems  are  non-stiff,  an  explicit/ 
implicit -predictor/corrector  method  may  be  employed  [28].'  Pre¬ 
dictor/corrector  algorithms  require  fewer  function  evaluations  than 
Hunga-Kutta  methods  and  allow  interpolation  between  solutions  with  no 
additional  function  evaluations.  This  is  important,  since  the  com¬ 
putational  overhead  in  evaluating  the  reduced  system  of  differential 
equations  is  significant.  The  explicit/implicit -predictor/corrector 
algorithm  DE/STEP,  INTRP  [  28  ]  is  used,  since  it  is  well  developed 
and  easy  to  adapt  to  the  equation  requirements.  In  addition,  the 
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integration  algorithm  DE  estimates  relative  stiffness  and  terminates 
execution  if  the  system  becomes  too  stiff. 

A  numerical  integration  algorithm  based  on  the  development  in 
previous  sections  and  the  explicit  integration  method  is  now  pre¬ 
sented.  The  following  steps  outline  the  numerical  integration 
algorithm . 

Let  i,  initially  0,  be  an  indicator  for  the  current  time  step; 
i.e.,  i  =  0  implies  t  =  t  .  Steps  1  to  5  are  satisfied  when  start¬ 
ing,  from  results  of  the  preceding  algorithm,  so  start  with  step  6. 

Step  1.  Check  j|$||  <  e.  If  satisfied,  the  extrapolated 
dependent  variables  u  are  satisfactory,  therefore  go  to  step  4. 

A  polynomial  extrapolator  is  maintained  for  u  by  the  integra¬ 
tion  algorithm  at  little  additional  overhead.  The  term  extrapolator 
is  used  rather  than  predictor,  because  the  program  does  not  maintain 
full  error  control  on  the  variables  u.  Since  this  extrapolator  may 
be  of  any  order  up  to  12,  the  constraints  $(u,  v,  t)  will  usually 
satisfy  the  above  closure  test,  thus  avoiding  the  following  steps  2 
and  3.  In  addition  the  partitioning  of  q  into  sets  u  and  v  depends 
upon  q  and  t  because  $  is  a  function  of  q  and  t.  Since  the  extra- 
polator  for  u  and  the  predictor  for  v  have  the  same  form,  it  is  not 
necessary  to  interrupt  the  integration  process  to  account  for  dif¬ 
ferent  sets  u  and  v.  Error  control  within  the  integration  algorithm 
is  maintained  only  on  v  since  steps  1  to  3  are  equivalent  to  main¬ 
taining  error  control  on  u. 
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Step  2.  Evaluate  and  perform  L-U  factorization  to  determine 
$  ,  5y,  H  and  the  partitioning  of  q  into  u  and  v. 

Step  3-  Iterate  to  determine  u,  using  Equation  3.5,  until 
11*11  <  «• 

Step  4.  Evaluate  8  ,  factor,  and  calculate  u  from  Equations 
3.16  and  3.17* 

The  partitioning  of  q  inco  u  and  s  s  v  depends  upon  q  and  t, 
because  is  a  function  of  q  and  t.  Since  all  elements  of  q  are 
candidates  for  s,  a  polynomial  extrapolator  for  u  and  a  predicator 
for  s  are  maintained  by  the  numerical  integration  algorithm,  with 
error  control  only  on  s.  This  facilitates  smooth  transition  of 
elements  of  q  into  and  out  of  s  within  the  integration  algorithm, 
with  no  interruption  of  the  integration  process.  Error  control 
within  the  integration  algorithm  is  maintained  only  on  s  since 
Equations  3.16  and  3.17  are  equivalent  to  maintaining  error  control 
on  u. 

Step  5.  Calculate  s  from  Equation  4.11. 

If  required,  calculate  ii  from  Equations  3.22  and  3.23,  and  h 

.  T  T  T 

from  Equation  3.30,  where  q  =  [u  ,  s  ]  and  v  =  s.  These  calcula¬ 
tions  are  not  required  for  the  integration  steps  that  follow,  however 
U  may  be  used  to  improve  the  extrapolator  of  u  [28]. 

Step  6.  Using  the  numerical  integration  PECE  algorithm  [28], 

predict  q  and  q  at  the  (i  +  l)St  time  step.  That  is,  predict  v^+1 

,  i+1  ,  .  ,  .  i+1  ,  .i+1 

and  s  and  extrapolate  u  and  u 

Step  7.  Execute  steps  1  to  5  to  evaluate  u,  v,  s,  and  U. 
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Step  8.  Using  u,  v,  s,  and  U  from  step  1 ,  correct  u,  v,  s, 

and  u. 

Correcting  u  and  u  is  done  only  to  improve  the  extrapolators 
of  u  and  u  for  the  next  step  which  reduces  predictor  error  if 
elements  of  u  and  v  are  subsequently  picked  up  by  v  and  s  at  the  next 
step. 

Step  9.  Estimate  integration  error  on  v  and  s.  Adjust  the 
current  time  step  and  integration  order  to  suit  integration  error 
requirements.  If  integration  error  is  too  large,  go  back  to  step  6, 
otherwise,  step  10. 

Step  10.  Execute  steps  1  to  5  again  to  obtain  updated  values 
for  q,  q,  and  k  and  report  the  solution  if  at  or  past  the  desired 
reporting  times.  Increment  i  and  return  to  step  6  or  stop,  if  the 
final  time  is  reached. 

4.5  Sparse  Matrix  Considerations 
(a)  Initial  Conditions 

As  noted  in  Chapter  2  the  nodal  method  of  constraint  formula¬ 
tion  yields  loosely  coupled  equations,  with  a  corresponding  sparse 
Jacobian  matrix.  If  the  system  has  a  large  number  of  constraints, 
sparse  matrix  manipulation  algorithms  will  be  efficient.  The  various 
equations  in  Chapters  3  and  4  can  be  expressed  in  different  forms  to 
avoid  unnecessary  calculati'--’  of  sparse  matrix  products.  The  pro¬ 
cedure  for  establishing  initial  conditions,  similar  to  Section  4.3, 
is  described  in  the  following  steps: 
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Step  1.  Specify  an  initial  estimate  q  of  the  constrained 
system  configuration.  This  estimate  must  be  near  the  expected  final 
configuration  to  insure  convergence.  If  static  equilibrium  is  to  be 
obtained,  the  estimate  must  be  near  the  desired  equilibrium  config¬ 
uration  to  avoid  convergence  to  alternate  stable  or  unstable 
configurations,  or  divergence.  Set  q  =  q^. 

Step  2.  If  some  elements  of  q  are  to  remain  fixed  during 
initial  assembly  of  the  constrained  system  or  static  equilibrium, 

append  to  Equation  4.4,  additional  algebraic  constraints  of  the  form 

0 

q .  -  q  .  =  0,  where  j  identifies  the  variables  to  be  held  fixed. 

3  3 

Other  types  of  algebraic  relations  between  coordinates  can  also  be 
formulated,  as  needed.  These  constraints  will  be  removed  prior  to 
step  6. 

Step  3.  Evaluate  $  (q)  (including  constraints  of  step  2)  and 

q 

factor,  as  described  in  Appendix  A,  to  determine  a  submatrix  of 
maximal  rank  and  to  determine  a  partitioning  of  q  into  u  and  v 
(dependent  and  independent  sets,  respectively).  If  does  not  have 
full  row  rank  then  too  many  constraints  or  inconsistent  constraints 
are  specified.  Therefore,  sane  equations  are  dependent  and  they 
can  be  temporarily  ignored  in  the  iterative  procedure  of  step  4.  A 
matrix  f  of  rank  3n  is  constructed  by  appending  the  matrix  I  to  }  , 

q  q 

as  done  in  Section  3. 5*  The  matrix  T  has  l's  in  columns  corresponding 
to  elements  v  of  q.  This  matrix  is  then  factored  using  sparse  matrix 
algorithms . 
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Step  4.  Iterate  using  Equation  3.3 6  to  solve  the  independent 
constraint  equations  of  step  2,  until  |  |f|  |  <  e,  where  €  is  a  spec¬ 
ified  closure  tolerance  level.  If  there  are  dependent  constraint 
equations,  check  them  for  violation.  When  dependent  constraints  are 
violated,  inconsistent  constraint  conditions  have  been  specified  and 
a  meaningless  solution  may  result.  If  static  equilibrium  is  not 
desired,  go  to  step  6. 

Step  5.  Determine  a  new  estimate  for  v  that  reduces  | |rV| |. 
From  Equation  3-39  expressed  as 

C[RV]  =  ‘  0 

evaluate  R  and  bv-  numerical  differencing  evaluate  dR  /dv.  Using 
the  method  of  Appendix  A,  factor  the  above  matrix  dRV/dv  and  solve 
the  equation 


dR 

dv 


Av  =  -  R 


for  Av.  If  | |RV| |  >  6,  a  specified  force  unbalance  tolerance,  return 
to  step  3.  If  the  iteration  in  steps  3  to  5  diverges,  either  incon¬ 
sistent  constraints  have  been  imposed  or  initial  estimates  of  system 
configuration  are  in  error. 

Step  6.  Specify  initial  estimates  for  velocities  q°  and  set 
q  =  q^.  The  velocities  q  are  adjusted  to  satisfy  Equation  3.13, 
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evaluated,  at  q  from  the  above  steps.  If  some  elements  of  q  are  to 
remain  fixed  at  q^,  append  to  Equation  3.13  additional  equations  of 
the  form  q  -  qP  =  0,  where  k  identifies  the  variables  held  fixed. 
Other  constraint  relations  between  velocities  can  be  formulated  as 
needed.  These  constraints  will  be  removed  prior  to  step  1  of  tran¬ 
sient  analysis. 

Step  7*  Factor  the  modified  matrix  of  Equation  3.13  as 
described  in  Appendix  A  to  determine  a  partitioning  of  q  into  u  and 
v  (not  necessarily  the  same  partitioning  as  in  step  3).  Append  to 
Equation  3*13  additional  equations  of  the  form  v  -  v°  =  0  to  form 

q  +  f  =  0.  This  is  equivalent  to  appending  a  matrix  T  to  § 
q  t  q 

with  l's  in  columns  corresponding  to  elements  v  of  q.  Factor  as 
in  step  3.  Calculate  u  using  the  above  partitioning  in  Equation  3.39* 
If  the  above  matrix  lacks  full  row  rank,  some  equations  are  dependent 
and  must  be  checked  for  consistency.  If  dependent  equations  are  not 
satisfied,  inconsistent  constraints  on  velocity  are  specified  and  a 
meaningless  solution  may  result.  Otherwise  a  consistent  set  of 
initial  conditions  is  now  available  for  transient  analysis. 

(b)  Numerical  Integration  of  the 
Equations  of  Motion 

The  steps  in  Section  4.4,  modified  to  incorporate  sparse 
matrix  techniques,  are  as  follows  (See  Section  4.4  for  additional 
comments) : 

Let  i,  initially  0,  be  an  indicator  for  the  current  time  step; 
i.e.,  i  =  0  implies  t  =  t^.  Steps  1  to  5  are  satisfied  when  starting 
from  results  of  the  preceding  algorithm. 
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Step  1.  Check  ||f)|  <  e.  If  satisfied,  the  extrapolated 

dependent  variables  u  are  satisfactory,  so  go  to  step  4. 

Step  2.  Periodically  evaluate  §  and  perform  L-U  factoriza- 

<1 

tion  to  determine  the  partitioning  of  q  into  u  and  v.  Append  to 
Equation  4.4  additional  equations  of  the  form  v  -  v^  =  0  to  form 
f  =  0  where  are  the  predicted  independent  variables.  Factor  the 
matrix  using  sparse  matrix  algorithms  and  go  to  step  3. 

Step  3.  Iterate  to  determine  u,  using  Equation  3.36,  until 

| |*| I  <  e. 

Step  4.  Evaluate  $  ,  factor  and  calculate  q  from  Equation 
•  T  T  T 

3.37*  Noting  that  q  =  [u  ,  s  ]  ,  this  step  is  equivalent  to  evalu¬ 
ating  the  differential  equations  v  =  s. 

Step  3.  Combine  Equations  2.19,  3.l8,  and  3.19  in  matrix 

form  as 


— 

“ 

“ 

M 

a  T 

*q 

q 

A  +  Q 

0 

-X 

-  C 

This  sparse  matrix  is  factored,  using  sparse  matrix  algorithms,  and 

the  accelerations  and  Lagrange  multipliers  are  determined.  Since 
..T  .T  T 

'4  =  [u  ,  s  ]  this  step  is  equivalent  to  evaluating  the  differential 
equations  for  s. 

Step  6.  Using  the  numerical  integration  PECE  algorithm,  pre¬ 
dict  q  and  q  at  the  (i  +  l)St  time  step.  That  is,  predict  v^+1 
si+1  and  extrapolate  u*+1  and  u**1. 


and 
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Step  7 .  Execute  steps  1  to  5  to  evaluate  u,  v,  s,  and  U. 

Step  8.  Using  u,  v,  s,  and  U  from  step  7 ,  correct  u,  v,  s, 

and  u. 

Step  9.  Estimate  integration  error  on  v  and  s.  Adjust  the 
current  time  step  and  integration  order  to  suit  integration  error 
requirements.  If  integration  error  is  too  large,  go  back  to  step  6, 
otherwise,  step  10. 

Step  10.  Execute  steps  1  to  5  again  to  obtain  updated  values 
for  q,  q,  and  X,  report  the  solution  if  at  or  past  the  desired 
reporting  times.  Increment  i  and  return  to  step  6  or  stop,  if  the 


final  time  is  reached. 
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CHAPTER  5 

PIECED  INTERVAL  ANALYSIS 

5.1  Introduction 

Dynamic  systems  characterized  by  one  or  more  events  of  short 
duration,  relative  to  total  system  response  time,  that  have  signifi¬ 
cant  effect  on  the  state  are  analyzed.  Apart  from  these  events,  the 
system  response  is  smooth.  Numerical  integration  algorithms  that 
descretize  the  time  domain  may  miss  events  such  as  impulsive  loads 
or  may  search  inefficiently  for  abrupt  state  transition  regions.  As 
transition  regions  approach  distinct  event  "logical  times",  discon¬ 
tinuities  occur  and  the  algorithms  will  fail  to  find  a  solution.  The 
alternative  is  to  sense  logical  times  and  handle  each  in  the  most 
efficient  and  accurate  manner  possible. 

Pieced  interval  analysis  implies  division  of  the  time  domain 
into  discrete  intervals  or  stages.  Stages  may  represent  discon¬ 
tinuous  logical  events  of  zero  time  duration,  such  as  momentum 
balance,  continuous  abrupt  transition  regions  of  short  time  duration, 
or  long  smooth  regions.  Each  stage  may  require  a  different  set  of 
state  equations,  initial  conditions,  and  analysis  technique. 

The  points  in  time  that  separate  stages  are  called  logical 
times.  They  are  determined  by  the  occurance  of  well-defined  events. 


81 


Since  these  logical  times  are  usually  functions  of  the  system  state, 
it  may  be  inpossible  to  specify  them  in  advance.  At  best,  one  may 
write  logical  functions  of  the  state  that  analytically  define  the 
occurance  of  important  events.  From  these  events  the  corresponding 
logical  times  are  obtained.  A  logical  events  monitor  built  into  the 
numerical  integration  algorithm  is  most  effective  in  locating  events 
and  obtaining  solutions  at  the  logical  times. 

Associated  with  the  logical  events  monitor,  hence  logical 
times,  is  user  selectable  programmed  control  logic  that  controls  the 
analysis  in  the  following  stage(s).  Actions  taken  depend  upon  the 
characteristics  of  each  given  problem  and  the  type  of  analysis  per¬ 
formed  within  each  stage.  Such  actions  may  be  as  simple  as  imposing 
limits  on  integration  time  stepsize  or  relaxing  error  tolerances 
during  short  duration  events,  such  as  impulsive  loading,  or  they  may 
be  as  extensive  as  obtaining  new  initial  conditions  through  momentum 
balance,  redefining  the  system  by  adding  or  deleting  bodies  and  con¬ 
straints,  selecting  the  most  appropriate  integration  algorithm,  and 
restarting  the  integration  process. 

An  important  application  of  pieced  interval  analysis  arrises 
when  dynamic  analysis  techniques  are  used  in  conjunction  with  design 
sensitivity  analysis  and  optimal  design  algorithms.  The  basic 
optimal  design  problem  with  fixed  starting  and  terminal  end  points, 
no  state  discontinuities,  and  no  intermediate  constraints  is  readily 
solved  by  existing  techniques  [  J>k],  Design  sensitivity  analysis 
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and  optimization  of  problems  with  discontinuities  at  logical  times 
requires  that  one  relax  these  restrictions  and  allow  for  variation 
in  logical  times  and  discontinuities  that  occur  there.  Such  methods 
have  been  used  for  small  scale  problems  [ 43  ] .  The  objective  here 
is  to  extend  and  apply  these  methods  to  large  scale  systems  with 
automated  equation  generation  and  solution. 

5.2  Logical  Events  Monitor  for 
Pieced  Interval  Analysis 

Logical  events,  defining  predictable  points  in  time  and  state, 
are  employed  to  facilitate  efficient  and  orderly  analysis  of  complex 
dynamic  mechanical  systems.  Boundaries  (marked  by  logical  events) 
of  successive  stages  in  analysis  often  must  be  located  before  they 
are  encountered  in  order  to  avoid  numerical  difficulties  associated 
with  the  prediction  of  system  state  variables  through  abrupt  or 
discontinuous  transition  regions .  Therefore  a  "logical  events 
monitor"  is  developed  that  utilizes  the  polynomial  predictor  of  the 
numerical  integration  algorithm,  to  solve  for  logical  times  before 
predicting  system  state  variables.  Having  identified  one  or  more 
logical  times,  the  monitor  forces  a  solution  precisely  at  the  first 
such  time,  interrupts  execution,  and  returns  control  to  the  user  or 
user  supplied  subroutines  for  further  action. 

To  define  the  boundaries  of  successive  stages  in  the  analysis 
of  a  mechanical  system,  a  set  of  logical  times  t^,  i  =  1,...,  k',  is 
defined  by  equations 
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=  fiJ'(q,  q,  t)  =  0,  j  =  1,...,  m'  (5 - 1) 

The  ordering  of  logical  times  is  defined  by  the  dynamical  system 
state  and  time. 

To  utilize  the  polynomial  predictor  of  the  numerical  integra¬ 
tion  algorithm.  Equation  5.1  is  differentiated  with  respect  to  time 

f  q  +  q  +  fij,  j  =  1,...,  m'  (5-2) 

q  q  t 

and  integrated  along  with  the  system  equations  of  motion.  This 
allows  freedom  in  the  selection  of  integration  algorithms,  because 
the  form  of  predictor  is  not  involved  in  the  formulation  of  Equation 
5.2.  Polynomial  predictors  use  past  history  of  state  variables  or 
their  derivatives  to  extrapolate  or  predict  them  ahead  in  time.  The 
logical  variables  i?  are  thus  predicted  ahead  in  time  before  the 
system  state  variables  are  predicted.  If  one  or  more  logical  vari¬ 
ables  should  pass  through  zero,  corresponding  logical  times  have  been 
passed.  Each  logical  time  is  then  found  by  solving  for  the  zero  of 
its  corresponding  polynomial  predictor.  The  event  corresponding  to 
the  smallest  time  step  is  thus  chosen  and  becomes  the  point  in  time 
at  which  a  solution  is  forced.  If  no  logical  variable  passes  through 
zero  in  a  given  time  step,  the  integration  algorithm  proceeds  as 
usual. 

The  success  of  this  method  relies  on  locating  the  zero  crossing 
of  each  logical  variable  i ^ ,  and  on  determining  a  reduced  stepsize 
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that  forces  a  solution  at  the  associated  logical  time.  It  is  expected 

that  logical  functions  are  smooth  enough  so  that  no  more  than  two 

zero  crossings  occur  within  any  time  step.  It  is  possible,  for 

example,  when  a  logical  function  defines  the  starting  and  ending 

times  of  a  very  short  impulse,  that  a  logical  function  passes  through 

zero  twice  within  a  single  time  step.  Each  zero  crossing  must  be 

located  in  succession  for  the  algorithms  to  work  properly. 

Presume  that  the  integration  algorithm  has  successfully 

th 

reached  time  station  t  and  the  j  logical  variable  has  the  value 

and  time  derivative  s^ .  The  logical  variable  t*  and  slope  s^  are 
n  n  n  n 

predicted  to  time  station  t^+^  as  and  s^+^.  Figure  5-1  illus¬ 

trates  various  possible  combinations  of  and  s^ .  First,  if 
i  1  11 

l  •  t  ,  >  0  and  s°  •  s°  .  >  0  there  are  no  zero  crossings.  Second, 
n  n+1  n  n+1 

if  t?  •  $  ,  >  0  and  s^  •  s^  <  0,  there  are  either  no,  one,  or  two 
n  n+1  n  n+1 

11  11 

zero  crossings.  Third,  if  j <  0  and  s^  •  s"+^  is  positive 
or  negative  there  is  one  zero  crossing.  Case  1  requires  no  further 
action.  Case  d  requires  interval  reduction  until  either  case  1  or 
case  3  is  achieved.  Case  3  requires  finding  the  zero  crossing  of 
and  taking  necessary  action  according  to  the  programmed  logic  assoc¬ 
iated  with  that  event. 

Finding  the  time  step  corresponding  to  the  zero  crossing  of 
j r  requires  finding  a  root  of  a  x  n  order  polynomial,  where  k  is  the 
order  of  the  predictor  employed  in  the  numerical  integration 
algorithm.  The  form  of  this  polynomial  depends  on  the  integration 
algorithm's  internal  representation  of  the  system  state  history.  Two 
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methods  sire  employed  in  many  predictor -corrector  integration  algor¬ 
ithms:  l)  a  Nordsieck  vector  representation,  which  stores  k  deriva¬ 
tives  of  the  state  variables  at  time  station  t  [  hk  ] ,  and  2)  a 
divided  difference  table,  based  on  the  t  and  previous  k  values  of 
the  state  veiriables  [  U5  ] .  Root  finding  techniques  for  each  method 
are  straight  forward  and  will  not  be  elaborated  upon  here. 

5-3  Intermittent  Motion  in  Constrained  Systems 

The  phrases  "discontinuous  forces",  "impulsive  loads",  and 
"discontinuous  motion"  appear  often  in  the  literature.  In  a  micro- 
mechanics  sense,  such  phenomena  do  not  exist.  However,  to  an 
observer  who  is  interested  in  the  dynamic  response  of  systems  over 
long  periods  of  time,  the  idea  of  "discontinuous"  or  "impulsive" 
may  be  quite  satisfactory.  On  the  other  hand,  an  observer  interested 
in  system  dynamic  response  over  short  periods  of  time,  specifically 
during  significant  events,  will  treat  them  in  a  continuous  manner. 
Discontinuous  and  impulsive  are  thus  labels  that  may  be  placed  on 
continuous  events,  which  in  a  macro -mechanics  sense  are  effectively 
represented  as  discontinuous  or  impulsive  in  the  observer's  time 
frame.  This  illustrates  the  desirability  of  pieced  interval  analysis 
for  certain  applications,  where  entirely  different  mathematical 
models  (and  possibly  analysis  techniques)  may  be  employed  during 
different  intervals  of  system  motion. 

In  this  section  momentum  balance  equations  are  developed  from 
Lagrange's  equations  of  motion,  assuming  that  discontinuous  external 
forces  or  impulses  are  applied  to  a  constrained  system.  The  equations 
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are  initially  written  for  the  case  in  which  these  events  occur  in  a 
finite  amount  of  time.  It  is  assumed  that  the  time  interval  ,  t^) 
is  small  enough  that  the  system  configuration  does  not  change 
appreciably,  that  is,  q(T2)  «  (qCT^).  In  addition,  the  constraints 
in  Equation  3.1  are  continuous  functions  of  q  and  t. 


T 


1 


t 


Figure  5*2.  Event  interval 


Let  t^  be  a  point  in  time  at  which  a  "violent  event"  occurs, 
which  is  to  be  approximated  by  a  discontinuity.  In  reality,  the 
event  occurs  over  a  time  interval  <  t  c  t^,  as  shown  in  Figure 
5*2,  and  behavior  is  smooth  except  possibly  at  t^.  Integrate 
Equation  4.8  to  obtain 


n  r.,W  uT,,UV  ,  /.,VU  „T, UUW1„ 

;  [M  +HM  +  (M  +HM  )Hjv  dt 
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-  /2  t*V  ♦  QV  *  HT(AU  a  «“) 

T1 

-  (M™  +  HTMUU)u*]dt  (5-3) 


Since  H  is  differentiable,  integration  by  parts  and  using  the 
mean  value  theorem  gives 

[MW  +  HTMUV  +  (MTO  +  HTMUU)H]v  |  2 

T1 

T2 

-  J  {  ~  [MW  +  HTMUV  +  (MW  +  HTMUU)H]v 
T1 


,  .v  irT.u  /.,vu  ,  „T,.uuv.*, 

+  A  +  H  A  -(M  +HM  )u  Jdt 

t  T2 

=  J  QV  dt  +  HT  J  QU  dt  (5.4) 


A 

Where  H  is  a  matrix  whose  elements  are  those  of  H  evaluated  at  points 
in  (tx,  t2). 

Treating  Q  as  impulsive  at  t^,  the  integrals  of  generalized 
force  are  "generalized  impulse", 


and 
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Taking  the  limit  in  Equation  5-^  as  -*  t^  and  "2  ->  t  ,  noting  from 

#  v  T  u  /  vu. 

Equations  3.19  and  3-22  that  li  is  bounded  so  A  +  H  A  -  (M 

+  H  M  )U  is  bounded,  yields  the  "impulse  momentum"  equation  at  t^ 

as 


[M 


yv 


,,T  uv 
H  M  + 


(M 


,vu 


HTMUU)H][v(t.+) 


v(t.-)] 


=  PV  + 


T  u 
HP 


(5-5) 


This  prescribes  the  velocity  jump  in  v  due  to  impulsively  applied 
loads . 


It  is  important  to  note  that  Equation  5*5  involves  impulse 
and  momentum  of  all  elements  of  the  mechanical  system.  This  is 
crucial,  since  the  bodies  making  up  the  system  interact  through  con¬ 
straints,  so  an  impulse -momentum  balance  relation  involving  only  the 
bodies  on  which  the  impulsive  force  acts  is  impossible.  Deriving  the 
relation  of  Equation  5*5  by  manual  calculation  would  be  extremely 
difficult  and  time  consuming.  One  of  the  strongest  points  of  the 
method  presented  here  is  the  automatic  assembly  of  the  coefficient 
matrices  of  Equation  5«5- 


Figure  5*3.  Impacting  bodies 


90 


For  impact  of  bodies  i  and  J,  as  shown  in  Figure  5.3,  a 
coefficient  of  restitution  e  provides  the  relative  velocity  relation 
in  direction  n  (where  n  is  taken  as  a  row  vector)  as 

°  L^V)  '  =  ‘  ‘  **<*,->]  <5.6) 

or  with  »Tc  R3”,  Equatior  5.6  may  be  „ltte„ 

^  *(t.+)  +  N“  =  -  e  [nv  v(ti-)  +  Nu  u(tr) 

or  using  Equation  3.17,  this  is 

NV  +  NU  H  v(t  +)  +  NU  u*(t,+) 

=  "  6  ^  +  nU  H  *(V>~e  nU  (5.7) 

Equation  3.16  implies  that  u*  is  a  function  only  of  q  and  t  and  con¬ 
straints  are  continuously  differentiable  in  q  and  t.  Therefore 
^  ('k^+)  -  u  (t^-)  and  Equation  5-7  reduces  to 

NV  +  NU  H  v(ti+)  =  -  e  NV  +  Nu  H1  v(t  ) 

L  J  i 

-  (e  +  1)NU  u*(t±-) 


(5.8) 
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The  generalized  impulse  of  the  force  f(t)  in  Figure  5.3  is 


P  =  J*  f(t)NTdt  =  pNT  (5-9) 

t.  -€ 

1 


where 


Ve 


p  =  J  f(t)dt 


V€ 


(5-10) 


T  T 

Defining  the  partitioning  P  =  pN  and  P  =  pNU  ,  Equation  5.3  gives 

[MVV  +  HTMUV  +  (M™  +  HTMUU)H][v(ti+)  -  ❖(ti-)] 

T  T 

=  p  [NV  +  HTNU  ]  (5.11) 

The  above  equations  are  put  into  matrix  form  for  solution  as 
follows.  Subtract  the  identity 


[nv  +  rAf]  v(ti-)  =  [nv  +  jAi]  v^-) 
from  Equation  5*8  and  define 
Av'i  =  vCt^)  -  v^-) 


(5.12) 
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The  matrix  equation  thus  becomes 

Mw  +  hTmuv  +  (Mvu  +  hTmuu)h 

(NV  +  NU  H) 


0 


-  +  l)[(NV  +  NU  H)v(ti-)  +  NU  u^(ti-)} 


(5.13) 


whose  solution  yields  the  desired  velocity  jumps  and  magnitude  of 
impulse  at  the  body  surfaces. 

5.h  Pieced  Interval  Computational  Algorithm 

The  development  in  this  and  previous  chapters  is  combined  into 

a  general  pieced  interval  analysis  algorithm.  In  the  discussion  to 

follow,  for  notational  convenience,  assume  that  logical  variable 

identifies  t^,  the  end  of  the  ith  stage.  Associated  with  is  a 

s  t 

set  of  system  state  equations  for  the  (i  +  l)  stage  (which  may  be 

dynamic  equations  of  motion  or  momentum  balance  equations).  In 

addition,  a  different  set  of  constraint  equations  (which  may 

include  boundary  constraints  .>n  velocity)  may  be  used  at  this  stage. 

1+1 

A  different  set  of  logical  functions  ft  ,  compatible  with  the  above 
constraint  and  state  equations,  may  also  be  introduced. 
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Prior  to  implementing  the  algorithm,  it  is  assumed  that  all 
possible  system  events  have  been  anticipated  and  appropriate  logical 
event  predictors  and  corresponding  constraint  and  state  equations 
have  been  formulated.  Knowledge  of  the  sequence  of  events  is  unneces¬ 
sary,  unless  this  knowledge  can  be  used  to  reduce  the  complexity  of 
equation  formulation.  Further,  it  is  assumed  that  no  event  takes 
place  prior  to  achieving  the  initial  static  equilibrium  configuration 
or  appropriate  initial  conditions  at  time  t  =  t  . 

The  computational  algorithm  is  as  follows: 

Step  1.  Obtain  initial  conditions  consistent  with  constraints, 
using  the  algorithm  of  Section  4.3. 

Step  2.  Select  the  appropriate  set  of  state  equations,  con¬ 
straints,  logical  functions,  etc.,  based  on  the  current  active  logical 
events  monitor  flag  or  continue  with  previous  equations. 

Step  3.  If  a  momentum  balance  is  required,  solve  for  velocity 
jumps  as  in  Section  5*3  and  return  to  Step  2. 

Step  4.  Using  the  techniques  described  in  Section  5-2  check 
for  active  logical  variables  (one  or  more  t?  passing  through  zero)  in 
the  next  time  step  to  be  attempted  by  the  numerical  integration 
algorithm.  Solve  for  time  steps  corresponding  to  the  zeros  of  poly¬ 
nomial  predictors  of  active  (if  any),  set  the  integration  time 
step  to  the  smallest  step  (first  occurring  l^)  and  set  the  corres¬ 
ponding  logical  events  monitor  flag.  Advance  the  solution  to  this 
point  in  time  by  executing  Step  5. 
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Step  5-  Perform  dynamic  analysis  using  the  algorithm  described 
in  Section  4.4.  if  the  solution  is  advanced  by  the  above  time  step 
return  to  Step  2,  unless  the  end  of  the  simulation  interval  is 
reached,  in  which  case  stop.  If  a  smaller  time  step  is  taken  to 
meet  integration  error  requirements,  reset  the  logical  events  moni¬ 
tor  flag  (since  the  event  has  not  yet  been  reached)  and  return  to 
Step  2. 


5.5  Sparse  Matrix  Considerations 
To  maintain  program  efficiency  when  mechanical  systems  become 
large,  sparse  matrix  algorithms  should  be  employed  and  matrix 
products  avoided.  The  methods  of  sparse  matrix  confutations  described 
in  previous  chapters  will  replace  the  corresponding  full  matrix 
operations  in  the  algorithm  of  Section  5.4.  It  is  possible  to  show 
that  Equation  5-13  can  be  expressed  in  a  form  suitable  for  sparse 
matrix  computations  by  first  writing  it  in  permuted  form  as 


— 

“ 

uu 

M 

uv 

M 

5  T 
*u 

uT 

N 

AUi 

vu 

M 

w 

M 

*  T 

V 

vT 

N 

AVi 

*u 

*v 

0 

0 

NU 

NV 

0 

0 

-P 

— 

_  _ 

95 


[  e  +  1]p':u  u^-)  +  NV  vC^-)} 


(5-lU) 


where  ^  is  a  vector  of  Lagrange  multipliers  .  This  equation  is 
equivale’:"  to  the  unpermuted  matrix  equation 


M  $ 


T 


$  0 

q 


N  0 


~ 

r 

A4± 

-  X 

-  p 

-{e  +  1}{N  q(t .  -)} 


(5.15) 


To  show  the  equivalence  between  Equation  5*13  and  5-1*+  carry  out 
the  matrix  products  in  Equation  5-1*+ 


MUU  Au^  +  MUV  Av^  -  X  -  NU  p  =  0 


(5.16) 


9 6 


VU  .  VV  .  T  ~  v1 

M  Ln±  +  M  Av^  -  9y  \  -  N  p  =  0  (5.17) 

5u  L\  +  ®v  L\  =  °  (5*l8) 

NU  Au±  +  NV  Av1=-{e  +  1}{NU  u(t±-)  +  NV  v^-)}  (5-19) 

Since 

•■Vv1  *  ‘vW’ +  mv*  ■ 0  < 

*u  *i  (ti+)  +  *vVV^  +  5t(ti+)  =  °  » 

and 


*t<V> 


W1 


then 


MW}  "  +  MW)  •  W)} =  0 

and  Equation  5*18  is  satisfied.  Equation  5.18  also  implies  that 

=  "  *u  1  Sv  L<ri  =  H  (5.20) 


Equation  5.19  is  simply  the  second  equation  of  Equation  5.15 

where  Equations  5.17  and  5.20  are  employed. 

, .  T  ~ 

Since  5^  is  nonsingular  in  Equation  5.16,  solve  for  \  and 

substitute  into  Equation  5.17,  thus 


^  Ain  +  Mw  Av±  +  $vT  5u"T(-  Muu  Ain  -  Av± 


T  T 

U  \  V 

+  N  p)  -  N  p  =  0 


Combining  terms  and  using  Equation  5.20  yields 


{MVV  +  HTMUV  +  (M^  +  HTMUU)H]  Av.  -  (NV  +  HTNU  >p  =  0 


which  is  the  first  equation  of  Equation  5.15.  Equation  5.15  is  thus 
solved  to  determine  changes  in  velocity  due  to  impact  between  bodies. 
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CHAPTER  6 
NUMERICAL  EXAMPLES 

6.1  Introduction 

Numerical  examples  are  presented  to  demonstrate  generality 
and  improved  program  efficiency,  resulting  from  the  analytical  methods 
developed  in  previous  chapters.  Generality  is  demonstrated  by  the 
ease  of  representing  models  of  complex  dynamic  mechanical  systems. 
Efficiency  is  demonstrated  by  contrasting  computer  simulation  times 
for  the  various  solution  methods. 

The  analytical  methods  developed  In  Chapters  5  and  4  demon¬ 
strate  the  greatest  improvement  in  program  efficiency  when  mechanical 
systems  are  heavily  constrained.  The  integration  algorithm  of 
Chapter  2  iteratively  solves  (6n  +  r)  differential  and  algebraic 
equations,  whereas  the  algorithm  of  Chapter  4  solves  6n  -  2r  dif¬ 
ferential  equations ;  n  being  the  number  of  system  rigid  bodies  and 
r  the  number  of  independent  algebraic  constraint  equations.  Defining 
system  degree  of  freedom  as 

dof  =  3n  -  r 

the  ratio  of  number  of  equations  for  the  two  methods  is 
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(6n  -  2r)/(6n  +  r)  =  (2  dof)/(9n  -  dof ) 

For  systems  with  many  bodies  and  a  few  degrees  of  freedom,  this  ratio 
is  small. 

The  three  degree -of -freedom  link  gear  multiplier  consisting 
of  twelve  bodies,  shown  in  Figure  6.1,  has  an  equation  ratio  of 
O.O57.  The  one  degree -of -freedom  Peaucellier  Lipkin  mechanism  in 
Figure  6.2  consists  of  8  bodies  and  has  an  equation  ratio  of  0.028. 

A  significant  reduction  in  matrix  and  numerical  integration  opera¬ 
tions  is  expected.  Increased  efficiency  also  results  from  the  use 
of  single  precision  computer  arithmetic,  reduced  number  of  corrector 
iterations,  and  larger  time  step  sizes. 

Operations  count  for  solution  of  sparse  matrix  problems  is 
of  the  order  N^'^,  where  Ii  is  the  matrix  dimension  [  46  ] .  The  ratio 
[(3n)/(9n  -  dof)]1’^  is  a  rough  estimate  of  the  reduction  in  opera¬ 
tions  count.  When  dof  «  3n  this  ratio  approaches  3^  =  0.17. 

As  dof  -*  3n  the  modified  constraint  matrix  (see  Equation  3.36) 
corresponding  to  the  numerator  of  this  ratio  approaches  an  identity 
matrix,  requiring  substantially  less  than  (3n)^*^  operations.  More 
than  a  10  to  1  reduction  in  simulation  time  for  the  above  mechanical 
system  simulations  was  observed. 

■As  the  equation  ratio  increases,  one  expects  less  gain  in 
program  efficiency.  For  example,  the  mechanical  system  described 
in  the  next  section,  consisting  of  12  bodies  with  24  degrees  of 
freedom,  has  an  equation  ratio  of  .596  and  a  corresponding  4  to  1 
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reduction  in  execution  time.  Note  that  the  above  equation  ratio  for 
both  methods  includes  5  additional  differential  equations  that 
describe  control  system  dynamics.  Another  mechanical  system  (not 
presented  in  this  report)  consisting  of  8  bodies  and  13  degrees 
of  freedom,  has  an  equation  ratio  of  .44  and  a  corresponding  7  to  1 
reduction  in  execution  time  [21].  Matrix  operation  counts  also 
play  important  roles  in  reduction  of  computer  time. 

6.2  Two  Articulated  Army  M-115 
Armored  Personnel  Carriers 

(a)  Vehicle  Pitch  Position  and 
* orce  Feedback  Control 

A  research  test  vehicle  shown  in  Figure  6.3  was  built  by  the 
Army  for  the  purpose  of  investigating  off -road  performance  of 
articulated  tracked  vehicles  [47,  48],  This  vehicle  is  selected 
for  simulation  here,  because  most  of  its  parameters  are  known,  some 
field  test  data  is  available  for  comparison  purposes,  and  its  pitch 
articulation  (while  negotiating  terrain  and  obstacles)  can  be 
adequately  represented  by  a  planar  model. 

The  test  vehicle  consists  of  two  identical  M-113  Army  armored 
personnel  carriers,  coupled  together  by  a  spherical  ball  socket  that 
allows  relative  vehicle  pitch,  yaw,  and  roll.  This  connection  is 
located  close  to  the  rear  of  the  front  vehicle,  to  give  the  rear 
vehicle  an  increased  moment  arm,  thus  allowing  it  to  remain  nearly 
horizontal  as  the  front  vehicle  is  elevated.  This  allows  the 
articulated  vehicle  to  negotiate  step  obstacles  up  to  5  feet  in  height. 


Relative  pitch  and  yaw  motion  of  the  system  in  control  ]  c-'i  by 
two  hydraulic  cylinders  acting  between  the  vehicle?  as  shown  in 
Figure  6.3.  Yaw  motion  occurs  when  one  cylinder  extend.'  and  tlv 
other  simultaneously  contracts.  Yaw  will  nol  be  considered  in  this 
investigation.  Therefore,  equal  actuator  lengths,  thus  equal  u- 
ator  exte-sion  and  contraction  rates,  will  be  assumed.  The  hydraulic, 
control  system  is  capable  of  maintaining  equal  actuator  extensic  :< 
rates,  while  negotiating  various  obstacles. 

Relative  pitch  and  yaw  can  be  controlled  by  the  vehicle 
operator,  via  a  control  stick.  If  the  stick  is  pushed  to  the  right 
or  left,  the  front  vehicle  yaws  to  the  right  or  left,  respectively. 

If  the  stick  is  pulled  back  or  pushed  forward,  the  front  vehicle 
pitches  up  or  down,  respectively.  Combinations  of  stick  movements 
yield  general  yaw  and  pitch  motions.  It  is  assumed  that  only  pitch 
command  signals  are  given  in  the  following  planar  simulations. 

For  the  operator  to  relate  control  stick  displacement  to 
relative  vehicular  pitch  displacement,  a  proportional  control  is 
required.  That  is,  pulling  the  stick  halfway  back  causes  the  vehicle 
to  pitch  up  to  5C$  of  its  maximum  allowed  displacement,  etc.  Pro¬ 
portional  control  is  provided  through  positional  feedback  to  the 
control  system  by  monitoring  the  hydraulic  actuator  extension  or 
contraction.  Proportional  yaw  control  is  obtained  by  a  similar 
process . 

Typical  operator  actions  for  negotiating  a  terrain  obstacle 
might  be  as  follows:  To  negotiate  an  obstacle  pull  the  stick  back 
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far  enough  such  that  the  front  unit  can  mount  the  obstacle.  Gradu¬ 
ally  drive  forward  while  adjusting  stick  position  so  that  the  vehicle 
maintains  maximum  conformity  to  the  terrain  and  obstacle.  This 
technique  generally  gives  the  smoothest  ride  and  better  vehicle 
traction  and  control.  However,  the  operator  must  continually  observe 
the  vehicle's  position  relative  to  terrain,  in  order  to  make  these 
adjustments.  This  requirement  may  limit  his  ability  to  perform  other 
tasks. 

Maximum  terrain  conformity  is  achieved  if  the  hydraulic 
actuators  are  removed  from  the  vehicles.  Since  this  isn't  possible, 
the  alternative  is  to  maintain  zero  pressure  within  the  actuators. 
This  is  accomplished  by  providing  force  feedback  from  the  hydraulic 
actuators  to  the  operator,  via  the  control  stick. 

Suppose  the  operator  pulls  back  on  the  control  stick  and  the 
vehicle  begins  to  pitch  up.  Pressure  builds  up  in  the  rod  end  of 
the  cylinder  and  a  voltage  from  a  pressure  transducer  is  generated. 
This  signal  is  amplified  and  applied  to  an  electro-mechanical  device 
that  pulls  forward  on  the  control  stick,  thus  increasing  its  resist¬ 
ance  to  the  operator's  applied  force.  This  gives  him  a  feel  for  the 
pressure  within  the  hydraulic  actuator.  He  has  simply  to  move  the 
control  stick  in  the  direction  it  wants  to  go.  In  fact,  he  may- 
remove  his  hand  and  the  control  stick  will  automatically  position 
itself  for  minimum  pressure  and  best  terrain  conformity.  Therefore, 
force  feedback  allows  the  operator  to  devote  more  of  his  attention 


to  other  duties. 
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The  control  system  cannot  respond  to  system  pressure  changes 
instantaneously,  because  of  inertia  and  hydraulic  system  compliance. 
Thus,  at  best,  it  can  only  maintain  a  zero  average  pressure  over  some 
time  interval  when  the  vehicles  experience  rapid  relative  pitch 
displacements.  In  fact,  pressure  transients  may  become  quite  high, 
as  evidenced  in  results  of  simulations.  A  complete  description  of 
the  electro-mechanical -hydraulic  control  system  model  is  given, 
following  a  description  of  the  vehicular  system  model. 

(b)  Elements  of  the  Articulated 
Vehicle  Mechanical  Model 

The  primary  components  of  an  M113  APC  are  the  chassis  and  ten 
roadwheels,  five  per  side.  The  roadwheels  are  connected  to  the 
chassis  by  roadarms  that  are  rigidly  attached  to  torsion  bars  that 
provide  chassis  suspension.  Each  vehicle  has  two  front  drive  and 
two  rear  idler  sprockets.  Segmented  tracks,  as  shown  in  Figure  6.3, 
pass  beneath  the  roadwheels  and  around  the  sprockets  to  form  contin¬ 
uous  loops.  Power  is  transferred  from  the  engine  to  drive  sprocket 
and  then  to  track  that  propels  the  vehicle.  The  roadwheels  simply 
roll  along  the  inner  track  surfaces.  Shock  absorbers  provide  damping 
in  the  front  and  rear  suspensions. 

A  representative  planar  model  is  established  from  the  above 
description.  It  is  cost  prohibitive  to  establish  a  model  that 
minutely  simulates  the  actual  system.  Therefore,  various  simpli¬ 
fications  can  be  made  where  minor  effects  are  ignored.  Since  motion 
in  a  plane  is  considered,  model  symmetry  is  employed  to  reduce 


problem  size.  It  is  supposed  that  opposing  roadwheels  experience 
equal  displacements,  hence  can  be  combined  into  single  bodies  having 
double  mass  and  rotational  inertia.  Suspension  stiffness  and  damping, 
track  properties,  etc.,  are  doubled. 

Track  dynamics  is  not  being  studied  here,  so  the  mass  and 
rotational  inertia  of  the  track  segments  beneath  roadwheels  is  equally 
distributed  among  them.  That  part  of  the  track  supported  by  the 
chassis  is  lumped  with  the  chassis  mass  and  moment  of  inertia.  Track 
supportive  effects  on  chassis  and  roadwheels  are  significant  and  will 
be  included  in  the  model. 

An  outline  drawing  of  the  model  is  shown  in  Figure  6.4.  This 
model  consists  of  12  bodies,  two  chassis  (Bll  and  B12)  and  ten  road- 
wheel  pairs  (31  to  BIO).  Roadarms,  represented  by  massless  links 
described  in  Chapter  2,  connect  roadwheels  to  the  chassis.  To 
simplify  a  track  model,  roadwheels  are  free  to  rotate,  but  only  small 
rotations  about  initial  angular  positions  are  considered.  Large 
rotational  displacement  of  roadwheels  is  prevented  by  the  track  model. 

Vertical  force  versus  vertical  roadwheel -displacement -suspen¬ 
sion  curves  and  vertical  force  versus  vertical  roadwheel -velocity¬ 
damping  curves  were  determined  experimentally  (Figures  6.5  and  6.6). 
These  forces  are  incorporated  into  the  model  employing  features  of 
the  standard  spring-damper -actuator  elements  described  in  Chapter  2. 
Elements  SI  to  S10,  connected  between  chassis  and  roadwheel  centers, 
are  given  arbitrarily  long  lengths;  i.e.,  1000  inches  so  that  their 
forces  relative  to  chassis  remain  essentially  vertical  for  all 
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roadwheel  positions.  Vertical  roadwheel  displacement  relative  to 

chassis  is  given  by  -  l  . ),  where  l  is  the  current  element 

ij  iJ 

length  and  £q  is  a  constant  reference  element  length.  Vertical 

id 

roadwheel  velocity,  relative  to  the  chassis,  is  given  by  (-  4,^). 
These  standard  variables  are  available  to  user  supplied  subroutines 


in  the  program.  Taking  k.  _  and  c  as  zero  in  Equation  2.9  yields 

J 

an  actuator  element  force  of 


Suspension  and  damping  force  functions  F  (jt0  -  H.  .)  and  Fd(-  SL 

S  ij 

are  supplied  by  cubic  spline  curve  fits  f  i+9 »  50]  to  the  data 
(straight  line  approximations  have  been  successfully  employed  as 
well)  [  51  ] .  The  actuator  force,  evaluated  as 


ij 


) 


=  F.Ur 


'id 


id 


£.  .) 

id 


+  Fd(- 


V 


is  automatically  added  to  the  system  generalized  forces. 

(c)  Elements  of  the  Track  and 
Terrain  Interaction  Model 

Track  forces  on  the  chassis  and  roadwheels  are  included 
through  spring-damper  elements  Sll  to  S22  (Figure  6.4).  It  is  assumed 
that  roadwheels  remain  in  contact  with  the  track  throughout  the  simu¬ 
lation.  Therefore,  elements  are  connected  directly  to  the  bottom  of 
the  roadwheels.  Elements  Sll,  Sl6,  S17,  and  S22  are  also  connected 
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directly  to  the  chassis  at  the  approximate  points  of  track  tangency 

to  the  drive  and  idler  sprockets.  These  elements  provide  track 

interaction  forces  between  the  chassis  and  roadwheels.  Each  element 

segment  has  an  effective  stiffness  determined  by  its  approximate 

length  and  a  track  stiffness  per  unit  length.  Internal  track  damping 

is  included  to  account  for  rubber  bushing,  rubber  track  pad,  and 

ground  dissipation.  These  coefficients  depend  nonlinear ly  on  track 

displacement,  since  track  segments  cannot  support  compression. 

Tracks  are  given  a  pretension  so  they  will  remain  tight  in 

various  dynamic  situations.  Each  element  Sll  to  S22  is  thus  given 

an  initial  pretension  of  10000  pounds.  In  order  to  determine  changes 

in  track  tension  (from  equilibrium  preload  levels),  due  to  dynamic 

loading,  it  is  necessary  to  have  reference  track  lengths  i  =  i 

ij 

at  equilibrium  for  each  of  the  above  elements.  These  lengths  cannot 

be  determined  when  the  model  is  established,  because  the  equilibrium 

configuration  is  unknown.  They  can  be  evaluated  directly  by  the 

program  as  equilibrium  is  established  as  follows:  Set  k.  .  =  c.  .  =  0 

ij  10 

and  Fq  =  10000  lb  in  Equation  2.9  for  each  element  Sll  to  S22. 
ij  1 

Constrain  angular  rotation  of  each  roadwheel,  using  initial  condition 

flags  on  input  data  cards  to  keep  element  attachment  points  at  the 

bottom  of  the  wheels.  Solve  for  equilibrium  ancl  set  each  element 

£q  to  l.  .,  which  then  becomes  the  reference  length  for  dynamic 
ij 

analysis.  Equation  2.9  for  the  dynamic  track  model  then  becomes 


hj  ~  [kij(*ij)(£ij  "  l0±5)  +  Cij(iij^ij  + 


where  k.  .(l.  .)  and  c. . (j i.  . )  are  determined  to  prevent  the  coefficient 
id'  XJ  1J  1J 

of  / l  .]  from  becoming  negative. 

Sij  1,5 

Roadwheels  must  be  prevented  from  penetrating  the  terrain 
surface.  That  is,  penetration  should  be  no  greater  than  that  achieved 
by  local  wheel,  track,  and  ground  deformation.  Support  generally 
occurs  at  one  or  more  points  on  the  lower  arc  of  a  wheel,  due  to 
direct  terrain  contact  or  indirect  track  support .  Support  reaction 
forces  may  not  act  vertically.  The  general  model  for  arbitrary 
terrain  is  rather  complex.  For  the  current  discussion,  assume  a 
level  terrain  (a  more  representative  model  for  irregular  terrain  is 
presented  later),  hence  roadwheel  support  reaction  forces  remain 
vertical.  These  forces  may  then  be  moved  to  the  wheel  centers  with¬ 
out  changing  the  model.  Since  roadwheels  are  circular,  the  terrain 
reference  elevation  can  be  shifted  up  to  the  wheel  centers  for  deter¬ 
mination  of  relative  wheel-terrain  elevation.  Global  location  of  the 
ith  roadwheel  center  is  given  by  coordinates  and  y_^  (see  Chapter  2). 
If  g.  is  global  terrain  elevation  and  r  is  roadwheel  radius,  then 
wheel  penetration  into  level  terrain  is 

Pi  =  (g±  +  rw  -  y.) 

If  wheel  reaction  force  is  developed  according  to  the  curve  in 

Figure  6.7,  the  force  F^p^  is  then  added  to  Q^1  (see  Equation  2.10). 

If  p^  s  0,  Fw(p^)  =  0  and  the  i^*1  wheel  is  not  touching  the  terrain. 

A  damping  force  D  (p  ,  p  ),  with  p.  =  g.  -  y, ,  is  also  added  to 
w  i  i  ill 


Qy1  to  include  roadwheel  rubber-,  track  rubber-,  and  ground -damping 

effects.  Damping  force  reduces  to  zero  as  p^  decreases  to  zero. 

When  vehicle  pitch  displacements  are  large,  other  points  on 

the  chassis  may  contact  ground.  The  most  likely  places  are  the  drive 

and  idler  sprockets,  labeled  A  to  D  in  Figure  6.4.  Sprocket  centers 

are  located  by  global  points  x  and  y  .  With  sprocket  radii  r  , 

pj  p3  s3 


one  may  obtain  support  and  damping  forces  as  above, 
are  then  included  in  generalized  forces  and  Q, 
bodies. 


These  forces 
for  the  proper 


As  noted  in  earlier  discussion,  the  usual  procedure  for 
representing  wheel-terrain  interference  is  to  monitor  the  distance 
between  the  bottom  point  on  the  wheel  and  the  point  on  the  terrain 
directly  below  it.  If,  as  in  Figure  6.8,  the  terrain  slopes  or  is 
irregular,  other  points  on  the  wheel  may  make  first  contact.  Thus, 
the  wheel  will  be  at  an  incorrect  height  and  the  support  force  on  the 
wheel  will  act  incorrectly.  A  modified  wheel  trajectory  may  be 
■  'obtained  by  moving  a  rigid  wheel  of  a  given  radius,  say  that  of  the 
roadwheel,  along  the  terrain  surface  to  obtain  a  trajectory  of  the 
lowest  point  on  the  wheel  [  52  j .  The  wheel  trajectory  for  the  terrain 
in  Figure  6.8  is  shown  in  Figure  6.9«  As  the  lowest  point  on  a 
wheel  follows  this  trajectory,  there  will  always  be  some  point  on  the 
wheel  just  touching  the  terrain,  but  there  will  be  no  terrain  penetra¬ 
tion.  This  is  illustrated  in  Figure  6.10.  For  the  purpose  of  this 
simulation,  it  is  sufficient  to  assume  that  the  resultant  support 
force  vector  acts  along  a  line  passing  through  the  point  of 
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wheel-terrain  contact  and  wheel  center.  Thus,  there  will  be  a  com¬ 
ponent  of  force  tending  to  impede  leftward  motion  of  the  wheel  and, 
if  not  otherwise  supported,  the  wheel  will  move  down  the  slope.  The 
angle  that  this  force  makes  with  respect  to  the  positive  y  axis  is 
shown  in  Figure  6.11.  The  curves  of  Figures  6.9  and  6.11  are  pro¬ 
vided  to  the  program  in  the  form  of  cubic  spline  functions. 

The  normal  and  frictional  components  of  force  acting  on  the 
wheel  remain  to  be  determined.  As  noted  earlier,  if  global 
coordinates  and  locate  the  center  of  a  roadwheel  and  if  g(x^) 
is  the  trajectory  height  beneath  the  wheel,  as  in  Figure  6.9,  then 
wheel  penetration  may  be  approximated  by 

Pi  =  Cg(xi)  +  rw  ‘  yiJ  (6.1) 

where  r^  is  the  roadwheel  radius.  Normal  force  angle,  from  the 
function  h(x)  of  Figure  6.11,  is 


Qi  =  h^xi^ 


(6.2) 


Figure  6.7  illustrates  a  typical  nonlinear  force -displacement  rela¬ 
tion  between  a  roadwheel  and  ground.  Let  this  curve  be  represented 
by  the  function  f(xi),  so  the  total  vertical  force  acting  on  the 
wheel,  as  shown  in  Figure  6.12,  becomes 


F  =  F  +  F_  =  f  (p ) 
s  ny  fy  ^ 


(6.3) 


where  F  and  F_  are  the  vertical  components  of  the  normal,  F  2  0 
ny  fy  r  n 

and  frictional  F^  forces,  respectively.  The  subscript  i  has  been 
dropped  for  convenience.  In  terms  of  the  angle  a  and  the  sign  con¬ 
ventions  of  Figure  6.12,  the  following  expressions  hold: 


F  =  -  F  sin  a 
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Kitaire  6.12.  Terrain-wheel  contact  forces. 
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ground  speed.  Let  v.  be  the  actual  velocity  of  this  point  on  body  j 
at  any  instant  t,  then  the  actual  displacement  of  this  point  is 

t 

S .  =  S„.  +  f  v.dt  (6.Q) 

d  Oj  J0  d  ’ 

Track  tension  is  developed  according  to  error  in  the  desired 
vehicle  position  and  velocity  as 


T. 

3 


-  Ka<sd 


s  .) 

rd 


K  (v.  - 
v  3 


v  ) 
r 


(6.10) 


where  the  constants  K,  and  K  are  chosen  to  characterize  the  vehicle 

d  r 

power  train.  Observe  that  T.  may  be  positive  (vehicle  lags  behind 

desired  position)  or  negative  (vehicle  is  ahead  of  desired  position). 

If  T.  >0,  T.  is  applied  to  track  elements  Sl6  and  S22.  If  T.  <0, 

-  T .  is  applied  to  track  elements  Sll  and  S17.  The  maximum  tractive 
J 

effort  is  limited  to  the  maximum  force  developed  by  a  slipping  track. 

The  sequence  of  steps  to  evaluate  track  driving  forces  for  the 

front  vehicle  are  described.  If  T .  >0,  set  F/  =  T..  Find  F  for 

3  d  s 

sprocket  C,  as  shown  in  Figure  6.4,  using  Equation  6.3.  Assume  slip 
>  0  (Figure  6.12),  since  the  vehicle  is  being  propelled  to  the  left. 

If  Fg  =  0,  the  sprocket  is  not  in  contact  with  the  terrain.  Therefore, 
no  frictional  tractive  force  is  developed  beneath  the  sprocket,  so 
the  entire  force  F#  is  transferred  to  track  element  Sl6. 

If  sprocket  C  is  in  contact  with  ground,  Fg  >  0.  Therefore, 
determine  the  frictional  tractive  force  F^  acting  between  sprocket 
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and  ground,  using  Equation  6.8.  Apply  F^  to  the  sprocket  at  the 
point  of  sprocket -terrain  contact,  directed  along  the  tangent  to  the 
terrain  surface  at  this  point.  The  track  driving  force  F/  is  now 
reduced  by  and  is  applied  to  track  element  Sl6.  The  process  is 
then  repeated  in  a  similar  manner  for  each  successive  roadwheel, 
rear  to  front.  Calculate  F  for  roadwheel  (B5)  from  Equation  6.3  and 
F  from  Equation  6.8.  Apply  F^  to  roadwheel  (B5)  at  the  point  of 
roadwheel -terrain  contact,  directed  along  the  tangent  to  the  terrain 
surface.  Again  reduce  F'  by  F^,  and  apply  this  force  to  element  S15. 
Continue  this  process  until  F'  has  been  reduced  to  zero,  or  until 
sprocket  D  has  been  passed.  If  F'  becomes  zero,  the  vehicle  is  not 
traction  limited.  If  F'  >0,  the  track  is  slipping  and  the  desired 
reference  position,  in  Equation  6.10,  cannot  be  attained;  so  set 

S  =  S.  -  v  t  and  S  =  S.  the  attained  vehicle  displacement. 

0..  J  r  rj  j 

If  T .  <  0  in  Equation  6.10,  the  process  is  reversed.  That  is, 
drive  forces  are  applied  from  front  to  rear.  First  check  sprocket 
D,  roadwheel  iBl),  etc.  for  contact  and  apply  appropriate  forces. 

The  vehicle  is  bei.  g  driven  to  the  right  in  this  case  and  slip  is 
assumed  to  the  left  (Figure  6.12),  hence  it  is  negative. 

The  significant  feature  of  this  model  is  that  track  slippage 
at  each  wheel  is  allowed  if  there  is  enough  force  differential  in 
the  track  to  overcome  local  friction  jiF  .  The  remaining  road  wheels, 
once  F '  becomes  zero,  are  assumed  to  have  zero  frictional  forces 


acting.  This  will  not  be  true  for  the  actual  system,  especially  on 


rugged  terrain.  However,  these  forces  will  be  directed  randomly, 
some  forward  and  some  rearward,  due  primarily  to  terrain  irregularity, 
and  will  not  contribute  significantly  to  vehicle  dynamics. 

(d)  Elements  of  the  Articulation 
Control  Model 

The  hydraulic  actuator  is  modeled  by  using  the  actuator  of 

element  S23.  Actuator  force  (fa  ),  displacement  and  velocity 

23 

Vg,  are  standard  state  variables.  These  variables  are  coupled  to  the 
control  system  equations  to  complete  the  model. 

Figure  6.13  contains  the  block  diagram  representation  of  the 
electro -hydraulic  control  system  [47,  48]  and  depicts  the  inter¬ 
relationship  between  the  servo  control  system  and  the  dynamic 
mechanical  model,  where: 

Fr(s)  =  operator  supplied  fore  and  aft  force  to  the  control 
stick,  lbp 

d^  =  length  of  control  stick,  inches 

0R(s)  =  angular  displacement  of  control  stick,  radians, 

j  0p 1  <  .436  radians 

T  (s)  =  control  stick  limiting  torque,  in-lb^, 

tc  =  control  system  time  constant,  sec 

KRp  =  potentiometer  gain,  volts/radian 

F.(s)  =  reference  voltage  signal  from  control  stick 

potentiometer,  volts 
Kp  =  amplifier  gain,  volts/volt 

z(s)  =  actuator  displacement,  inches 
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=  total  actuator  length,  inches 
=  reference  actuator  length,  inches 
=  actuator  velocity,  in/ sec 

=  actuator  potentiometer  feedback  gain,  volts/in 
=  error  signal  voltage,  volts 
=  amplifier  gain  ma/volt 
=  amplifier  current,  ma 
=  servo  valve  flow  gradient,  GPM/ma 
=  pump  stroke  mechanism  velocity  gradient, 

%  stroke /sec 

=  servo  valve  time  constant,  sec 

=  pump  yoke  displacement,  %  stroke 

=  pump  yoke  limiting  force,  lb^/in 

=  pump  flow  gradient,  (i r?/sec)/{%  stroke) 

=  pump  flow,  irP/sec 

=  load  flow,  in^/sec 

=  leakage  flow,  in^/sec 

=  flow  due  to  compressibility,  in^/sec 

2 

=  actuator  area  -  rod  end,  in 

2 

=  actuator  area  -  cylinder  end,  in 

=  effective  bulk  modulus  including  system  compliance, 
psi 

=  entrained  volume  of  hydraulic  fluid 
=  combined  system  leakage  coefficient,  (in'’/sec)/psi 
=  system  pressure  at  the  actuator,  psi 
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=  force  developed  in  the  actuator,  lb^ 

=  pressure  transducer  gain,  mv/psi 
=  control  stick  force  feedback  gain,  in-lb^/mv. 

The  articulated  vehicle  control  system  may  be  divided  into 
three  major  blocks.  As  shown  in  Figure  6.13,  block  1  generates  a 
reference  voltage  that  ultimately  leads  to  a  given  hydraulic  actuator 
displacement.  The  major  component  in  this  block  is  the  mechanical 
hardware  that  interfaces  with  the  operator.  Block  two  contains  the 
hydraulic  pump  yoke  and  servo  valve  control  system.  In  response  to 
an  error  signal  from  block  1,  this  system  controls  direction  and 
fluid  flowrate  in  the  hydraulic  pump.  Block  3  characterizes  dynamic 
response  of  the  hydraulic  system,  including  the  actuator  connected 
to  the  articulated  vehicle  model. 

The  control  stick,  shown  in  Figure  6. it  [ !i7  ],  rotates  about 
a  fixed  point.  A  direct  current  linear  actuator,  shown  in  its  neutral 
position,  is  coupled  to  an  extension  of  the  control  stick  by 
springs  and  .  These  springs  provide  restoring  moments  to  main¬ 
tain  the  control  stick  in  its  neutral  position.  In  response  to 
pressure  buildup  in  the  hydraulic  actuator,  a  voltage  is  developed 
across  the  linear  actuator,  causing  springs  and  K to  shift  either 
left  or  right.  An  increased  moment  is  then  developed  on  the  control 
stick  to  oppose  the  operator  supplied  force,  thus  providing  force 
feedback  to  the  operator.  In  addition,  mechanical  stops  prevent 
more  than  +  25°  rotation  of  the  control  stick.  The  dynamics  of  this 
system,  including  inertial  effects,  is  included  in  block  1.  Control 
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stick  response  is  determined  by  three  moments:  operator's  moment, 

(Fd  '  d  ) ;  pressure  feedback  moment,  ( -K  K  P);  and  mechanical 
K  1  p  LA 

stop  limiting  torque,  [-  Tr(0_,)]. 

L  K 

A  rotary  potentiometer  connected  to  the  control  stick  develops 
a  reference  voltage  that  is  proportional  to  angular  rotation  0R. 

This  signal  is  then  amplified  and  combined  with  a  hydraulic  actuator 
displacement  signal  (-  z)  to  obtain  an  error  signal  E.  This 
error  signal  is  then  amplified  and  passed  to  the  pump  yoke  control 
system. 

Servo  valve  and  pump  yoke  characteristics  are  illustrated  in 
block  2  of  Figure  6.13.  Significant  nonlinear  effects  in  the  system 


3(t) 


forward 


Figure  6.l4.  Contro]  stick  and  force  actuator. 
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arc  saturation  of  the  servo  valve  flow  control  K  and  mechanical 

v 

limitations  on  pump  yoke  displacement  F  .  Inertial  effects  in  this 
block  introduce  additional  delay  in  overall  system  response. 

Block  three  outlines  the  hydraulic  flow  model.  Flow  rates 
within  the  system  are  moderate,  so  pressure  drop  within  the  lines 
is  small.  This  simplifies  the  model  significantly.  The  pump 


supplies  fluid  volume  flow  rate  Q  ,  to  make  up  for  leakage  losses 
G  ,  compressibility  and  actuator  volume  displacement  rate 

The  system  contains  a  significant  amount  of  flexible  tubing, 
which  results  in  reduced  system  compliance,  hence  a  small  effective 


bulk  modulus.  Fressure  limiting  is  obtained  by  further  reduction  of 
the  effective  bulk  modulus,  as  relief  valve  blowoff  begins.  Actuator 
flow  rate  from  the  pump  is  determined  by  the  product  of  actuator 
cross-sectional  area  A  and  velocity  v.  The  area  A. ^  on  the  rod  end 
of  the  actuator  is  smaller  than  the  piston  end  area  Ag.  Thus,  area 
is  adjusted  appropriately  as  a  function  of  actuator  velocity  to 
achieve  the  proper  flow  rate.  Force  developed  in  the  actuator  also 
depends  on  which  side  of  the  piston  pressure  is  developed.  Thus, 
area  is  adjusted  as  a  function  of  pressure,  for  the  purpose  of 
calculating  actuator  force.  The  controller  consists  of  two  identical 


hydraulic  systems,  each  driving  an  actuator.  In  this  simulation, 


only  one  hydraulic  system  is  modeled,  so  the  effective  actuator 


force  is  doubled  to  account  for  two  actuators. 
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(e)  Crossing  a  Three -Foot  Obstacle  - 
Numerical  Results 

The  three -foot  step  obstacle  shown  in  Figure  6.8  is  selected 
to  test  the  program's  ability  to  adequately  describe  system  dynamics. 
In  addition,  the  effect  of  force  and  position  feedback  on  vehicle 
performance,  while  negotiating  a  significant  obstacle,  is  to  be 
investigated.  A  velocity  of  40  inches/second  (2.27  miles/hour)  is 
selected  as  the  obstacle  crossing  speed.  This  is  close  to  the 
design  specifications  of  minimum  2.5  miles/hour  over  a  2.5  foot 
obstacle . 

Contrary  to  usual  practice,  the  vehicle  is  not  put  into  an 
initial  pitch  up  attitude  before  encountering  the  bank,  but  is  forced 
to  propel  itself  up  the  side  of  the  bank.  The  vehicle  is  run  in  a 
hands  off  mode,  because  it  is  difficult  to  program  operator  response 
to  the  complex  dynamic  environment. 

Three  simulations  of  obstacle  crossing  are  performed.  The 
first  case  is  with  force  and  position  feedback  active,  the  second  is 
with  force  feedback  deactivated,  and  the  third  is  with  the  entire 
control  system  deactivated.  In  the  last  situation  the  control  system 
is  passive,  except  that  the  hydraulic  pump  is  providing  enough  fluid 
flow  rate  to  replace  leakage  losses.  A  sequence  of  computer  generated 
line  drawings  of  vehicle  position  at  one -second  intervals  for  the 
three  cases  is  shown  in  Figures  6.15-6.17. 

The  following  observations  can  be  made  from  these  figures. 

Force  and  position  feedback  in  Figure  6.15  tend  to  minimize  actuator 


Figure  6.17.  Obstacle  crossing  with  the  control  system 
deactivated  -  ten  second  simulation  at  ijO 
inches/second. 
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stiffening  of  the  coupled  vehicle  system.  This  will  be  the  most 
comfortable  and  stable  situation  since  the  vehicles  tend  to  conform 
to  the  terrain. 

Position  feedback  alone  in  Figure  6.1 6  tends  to  maximize 
vehicle  rigidity  in  the  hands  off  mode  because  the  control  system 
continually  works  to  maintain  minimum  departure  from  the  zero 
relative  pitch  position.  With  operator  assistance,  however,  one 
might  expect  vehicle  response  approaching  that  of  case  1. 

The  passfv  model  in  Figure  6.17  is  slightly  more  flexible 
than  case  2,  because  of  system  compliance.  An  ideal  passive  system 
for  large  obstacle  crossing  would  be  one  that  allows  free  or 
partially  restricted  fluid  flow  in  the  actuator,  through  a  simple 
bypass . 

To  understand  more  clearly  how  force  feedback  works  to 
minimize  actuator  pressure,  it  is  shown  plotted  versus  time  in 
Figure  6.l8.  Note  that  positive  pressure  implies  that  pressure  acts 
on  the  piston  end  of  the  cylinder  and  is  zero  on  the  rod  end. 

Negative  pressure  implies  that  pressure  acts  on  the  rod  end  of  the 
cylinder  and  is  zero  on  the  piston  end.  This  is  a  convenient  way 
to  mathematically  represent  system  pressures.  Pressure  is  generally 
positive  for  the  first  3  to  4  seconds  (pressure  on  piston  end) ,  as 
the  front  vehicle  is  forced  up  and  over  the  step.  It  then  remains 
negative  between  4  and  10  seconds  for  case  2,  since  the  front  vehicle 
is  now  being  supported  by  the  actuator.  Pressure  is  essentially 
negative  between  4  and  8  seconds  and  positive  between  8  and  10 
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seconds,  when  force  feedback  is  active.  It  is  clear  that  force 
feedback  substantially  reduces  actuator  pressure,  as  is  expected. 


TIME  (SECONDS) 

Figure  6. 18.  System  hydraulic  actuator  pressure. 


The  reason  for  reduced  pressure  is  made  clear  by  looking  at 
Figure  6.19,  angular  displacement  of  the  control  stick,  and  Figure 
6.20,  hydraulic  pump  yoke  displacement.  When  actuator  pressure 
initially  goes  positive,  the  control  stick  is  driven  into  the  full 
pitch  up  position  (-  25°).  The  control-stick  signal  then  induces  a 
large  pump  yoke  displacement,  causing  hydraulic  fluid  flow  from  the 
high  pressure  end.  When  pressure  switches  to  the  opposite  side  of 
the  piston,  control  stick  and  pump  yoke  displacements  quickly  switch 
to  their  opposite  extremes,  again  to  compensate  for  large  actuator 


pressure . 


FORCE  MID  POSITION  FEEDBACK  ACTIVE 


Figure  6.19.  Angular  displacement  of  the  control  stick 
(negative  =  pitch -up ) . 
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(f)  Model  Verification  and  Discussion 

The  articulated  vehicular  system  simulated  in  this  study  was 
assembled  and  tested  by  the  Army  in  the  early  to  mid  1970's  [47, 

43  ]  .  A  large  amount  of  field  and  design  test  data  was  accumulated. 
Upon  termination  of  the  project,  the  vehicles  were  disassembled  and 
returned  to  regular  service.  Several  years  later  the  DADS  computer 
program  was  developed  and  it  was  recognized  that  this  vehicle  would 
make  an  excellent  test  model  to  demonstrate  DADS  analysis  and  design 
potential.  Unfortunately,  some  important  system  parameters  were  not 
measured  or  were  not  documented.  Therefore,  it  was  necessary  to 
estimate  some  parameters  in  order  to  complete  the  model. 

Several  tests  that  were  documented  [47,  48  j  measured  tran¬ 
sient  pitch  response  to  various  pitch-up  and  pitch-down  command 
signals.  The  test  vehicles  were  stationary  and  on  a  level  surface. 
Simulated  vehicle  pitch  rates  and  hydraulic  pressures  under  identical 
conditions  compared  favorably  with  field  test  data  [53]*  Discrep¬ 
ancies,  however,  did  provide  suggestions  for  improvement  in  the 
initial  model,  i.e.,  the  need  for  increased  hydraulic  system 
compliance,  track  supportive  effects,  and  frictional  effects.  It  is 
believed  that  these  initial  simulations,  for' which  field  test  data 
was  available,  have  resulted  in  a  more  representative  mathematical 
model. 

Unfortunately,  such  extensive  test  data  has  not  been  recorded 
for  obstacle  crossing  experiments.  However,  16  mm  filmed  records 
were  made  of  various  obstacle  crossing  tests.  One  such  test  involved 
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a  three  foot  step  obstacle  where  the  vehicle  was  placed  in  an  initial 
pitch -up  attitude  allowing  the  front  end  to  negotiate  the  obstacle. 
The  vehicle  then  propelled  itself  onto  the  step  with  some  operator 
assist.  Test  vehicle  velocity  was  comparable  to  the  velocity  in 
this  simulation. 

A  l6  mm  computer  generated,  animated,  film  of  the  dynamic 
simulation  (composed  of  images  similar  to  those  in  Figures  6.15- 
6.17)  was  compared  to  the  filmed  test  results  mentioned  above.  The 
agreement  was  significant.  Dynamic  response  of  the  vehicles  as  the 
roadwheels  encountered  the  step  was  very  similar.  The  time  lag  in 
the  vehicle's  ability  to  conform  to  terrain  was  apparent  in  both 
cases,  but  not  as  significant  in  the  field  test  results.  This  could 
be  due  to  operator  assist  or  to  less  inertia  in  the  actual  control 
system  than  was  modeled,  allowing  faster  system  response. 

This  example  is  presented  in  detail  to  emphasize  the  many 
important  considerations  one  must  make  in  order  to  have  even  a  simple 
representative  model  of  a  complex  system.  It  is  shown  that  most  of 
the  nonstandard  effects  are  conveniently  introduced  by  modifying 
existing  standard  program,  elements  and  employing  standard  state 
variables.  Thus,  the  magnitude  of  equation  development  and  pro¬ 
gramming  effort  has  been  kept  to  a  bare  minimum.  Even  more 
significant  than  the  efficiency  gained  by  the  analysis  methods 
developed  herein  is  the  fact  that  derivatives  of  generalized  forces 
are  not  required.  The  program  described  in  Chapter  2  requires 


139 


derivatives  of  generalized  forces,  which  substantially  increases  the 
amount  of  program  preparation  and  debugging  effort. 

The  independent  generalized  coordinates  identified  by  the 
program  for  the  above  simulations  are  x^,  x^_,  y^  to  y^’  to 

<p  Q.  One  could  make  intuitive  arguments  based  on  mechanical  or 
mathematical  principles  as  to  why  certain  variables  are  selected 
over  others.  However,  since  the  selection  process  depends  so  strongly 
on  constraints,  one  must  thoroughly  understand  each  constraint  type 
and  how  it  interacts  with  others,  if  an  acceptable  set  is  to  be 
selected.  This  emphasizes  the  importance  of  automatic  variable 
selection. 

It  is  emphasized  that  various  vehicle  and  control  system 
parameters  are  estimated,  because  measured  values  are  not  available. 
Thus  simulated  response  will  vary  somewhat  from  actual  system 
response.  The  primary  purpose  of  this  research  is  to  demonstrate 
that  DADS  can  be  effectively  employed  to  obtain  dynamic  response  of 
complex  mechanical  systems.  Therefore,  it  is  demonstrated  to  be 
an  invaluable  analysis  and  design  tool. 

6.3  A  Mechanism  With  Intermittent  Motion 
(a)  Introduction 

A  precision  weapon  mechanism  was  selected  for  simulation  to 
illustrate  use  of  the  pieced  interval  analysis  methods  developed  in 
Chapter  5*  Previously,  methods  of  intermittent  motion  analysis  have 
used  pieced  interval  analysis,  in  which  the  analyst  writes  the 
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equations  of  motion  between  times  at  which  discontinuous  events  occur 
[43] .  Momentum  balance  equations  were  also  written  to  account  for 
velocity  discontinuities  that  occur  in  various  system  configurations. 
Numerical  integration  was  halted  at  points  of  discontinuity,  new 
initial  conditions  were  formulated,  and  integration  was  then  con¬ 
tinued.  A  limitation  of  this  method  of  analysis  is  the  level  of 
effort  required  to  write  system  equations  that  are  valid  in  intervals 
between  events  whose  ordering  is  not  generally  known  before  the 
analysis  is  begun.  Thus,  the  analyst  is  required  to  write  equations 
and  computer  code  for  all  ordering  of  logical  events  that  may  con¬ 
ceivably  occur. 

One  method  that  has  been  used  to  alleviate  the  foregoing 
difficulty  is  to  use  Heaviside  step  functions  that  define  logic 
associated  with  the  events  occurring  during  intermittent  motion. 

These  discontinuous  functions  may  then  be  smoothed  to  provide  a  set 
of  governing  differential  equations  of  motion  [54],  This  pro¬ 
cedure  can  be  justified  on  the  basis  of  distribution  theory  [55, 

56]  and  has  been  successfully  employed  in  weapon  mechanism  dynamics 
[57].  The  distribution  theoretic  method  has  been  used  in  conjunc¬ 
tion  with  the  DADS  computer  program  to  automatically  generate  the 
system  equations  of  motion  [58],  by  defining  "logical  spring- 
dampers"  that  account  for  certain  aspects  of  intermittent  motion 
[21] .  These  analyses  are  plagued  by  integration  failure  and 
inefficiency,  while  integrating  through  abrupt  transition  regions. 
This  deficiency  has  now  been  eliminated  by  employing  the  logical 


events  predictor  developed  in  Chapter  5  to  locate  such  events  before 
they  are  encountered.  Thus,  the  integration  algorithm  can  be  pre¬ 
pared  to  handle  the  events  or  alternate  integration  methods  can  be 
employed. 

Alternately,  one  may  choose  to  replace  (or  supplement)  the 
smoothed  continuous  events  by  discontinuous  events,  in  -which  momentum 
balance  is  employed.  Costly  numerical  integration  through  transition 
regions  is  then  replaced  by  efficient  algebraic  solution  of  momentum 
balance  equations.  Following  momentum  balance,  one  may  elect  to 
retain  the  continuous  elements,  if  desired,  to  represent  interference 
between  adjacent  bodies. 

(b)  Description  of  an  Automatic  Cannon  System 

An  automatic  weapon  mechanism  shown  in  Figure  6.21  consists 
of  four  main  masses;  the  receiver  R  (assumed  rigidly  attached  to  the 
inertial  reference  frame),  the  barrel  assembly  B,  the  sleeve  S,  and 
the  sear  SR.  The  sequence  of  system  operation  for  a  cycle  is  as 
follows : 

The  barrel  assembly  is  initially  moving  to  the  right  at  40 
inches  per  second  at  the  instant  shown  in  Figure  6.21.  A  round  is 
inserted  into  the  shaded  region  of  the  chamber  as  the  barrel  is 
moving  forward.  A  link  pq  is  connected  at  point  q  to  a  sleeve  that 
slides  along  the  outer  surface  of  the  barrel  and  to  a  pin  p  located 
at  the  intersection  of  the  two  cam  slots  to  control  sleeve  position  on 
the  barrel.  The  receiver-cam  is  attached  to  the  receiver  and  the 
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barrel-cam  is  attached  to  the  barrel.  As  the  barrel  moves,  the  point 
of  cam  slot  intersection  changes,  forcing  the  pin  to  mov< • . 

Figure  6.22  shows  the  configuration  when  the  pin  just  begins 
to  move  upward  as  the  barrel  is  moving  to  the  right.  This  upward 
motion  drives  the  sleeve  to  the  left,  relative  to  the  barrel,  hence 
encasing  the  round.  Eventually  (Figure  6.23),  the  pin  reaches  the 
highest  point  on  the  receiver  cam,  and  the  sleeve  in  its  rearmost 
position  relative  to  the  barrel,  securely  locks  the  sleeve  in  posi¬ 
tion  to  form  the  chamber  of  the  weapon.  The  sleeve  now  moves  with 
the  barrel. 

Next,  as  shown  in  Figure  6.24,  the  barrel  contacts  a  front 
buffer.  This  buffer  is  designed  to  bring  the  barrel  safely  to  rest 
in  the  event  of  a  round  misfire.  Shortly  thereafter  (Figure  6.25) 
the  round  is  fired,  developing  an  impulse  on  the  barrel  that  reverses 
its  direction  of  travel. 

The  sequence  of  operations  is  then  reversed  (Figures  6.24, 
6.23,  and  6.22).  After  the  sleeve  has  been  retracted,  the  spent 
round  is  ejected  and  the  barrel  impacts  the  sear  (Figure  6.2l).  The 
sear  and  barrel,  moving  rearward  together,  are  then  brough  to  rest 
by  a  rear  buffer  (Figure  6.26).  The  sear  and  barrel  are  then  driven 
forward  by  the  buffer  spring  until  the  sear  impacts  a  stop  on  the 
receiver  (Figure  6.2l).  The  barrel  separates  from  the  sear  and  a 
cycle  is  complete  (ready  for  the  next  round  in  the  automatic  fire 
mode),  or  the  barrel  remains  attached  to  the  sear  to  terminate 


execution . 
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(c)  Elements  of  the  Model 

The  weapon  mechanism  consists  of  five  rigid  bodies:  receiver, 
barrel,  sear,  sleeve,  and  pin,  numbered  in  this  order.  The  sear  and 
barrel  each  move  on  the  receiver,  along  translational  joints,  and 
the  sleeve  moves  on  the  barrel  along  a  translational  joint.  Mass 
of  the  barrel,  sleeve,  and  sear  are  taken  as  2.919*  .2368,  and  . 2044- 
slug  -feet/inch,  respectively.  Other  masses  and  rotational  inertias 
are  ignored  or  are  not  important,  since  no  rotation  is  allowed.  The 
sleeve  is  connected  to  the  pin  by  a  massless  link  pq.  For  the  barrel 
position  shown  in  Figure  6.21,  the  link  is  initially  at  a  45°  angle 
and  it  moves  to  a  horizontal  position  in  Figure  6.23,  as  the  pin 
rides  up  the  cam  slots. 

The  pin-cam  constraints  are  represented  by  cubic  spline  curve 
fits  to  digitized  data  taken  from  the  mechanism.  The  vertical 
position  of  the  pin  relative  to  the  receiver  is  written  as  a  function 
of  the  horizontal  position  of  the  pin  relative  to  the  receiver  as 

*P  "  yR  _  fR(XF  -  XR}  =  °  (6-ll} 

where  f  is  the  functional  spline  representation  of  the  cam  path  on 

n 

the  receiver.  Likewise,  a  cam  path  on  the  barrel  is  written  as 


fB<*P  -  V  *  0 


(6.12) 
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These  equations  and  their  first  and  second  time  derivatives  are 
appended  to  the  standard  constraint  formulation,  to  complete  the 
model. 

Two  forces,  F^  and  F^,  drive  the  barrel  during  its  forward 

(counter  recoil)  and  rearward  (recoil)  motion,  respectively.  A  front 

buffer  B„  and  a  rear  buffer  B  slow  the  barrel  assembly  during 
f  r 

extreme  displacement.  Both  front  and  rear  buffers  are  designed  to 
produce  constant  retarding  forces. 

(d)  Logical  Events 

Logical  times  t^  at  which  impact  or  other  irregularities  of 
intermittent  motion  occur  are  introduced  as  an  integral  element  of 
the  dynamic  model.  Between  these  times,  the  motion  and  acceleration 
of  the  system  is  continuous.  At  these  times,  discontinuities  in 
velocities  and  acceleration,  changes  in  system  constraints,  and  mass 
capture  or  release  can  occur.  These  logical  times  are  functions  of 
the  system  state  and  are  determined  as  the  simulation  progresses. 
Logical  times  will  now  be  defined  for  the  firing  from  run-out  mode 
of  weapon  operation: 

(l)  tQ  =  0  (Figure  6.2l):  The  barrel  assembly  B,  in  auto¬ 
matic  fire,  passes  the  sear  position  with  velocity  x2  = 
40  in/sec.  Initial  starting  point  is  not  considered  as 
a  logical  event.  A  forward  driving  force  Ff  =  1600  lb^ 


acts  on  the  barrel. 


151 


(2)  t  (Figure  6.24):  The  barrel  B  contacts  the  front  butter 
and  =  -  69OO  lb^  becomes  active.  Restart  integration 
because  of  discontinuous  acceleration. 

(3)  tg  (Figure  6.25):  The  charge  is  ignited.  An  impulse  of 

-  880  lb  -sec  acts  on  the  barrel  B.  Perform  momentum 

balance  to  obtain  new  velocities.  F  is  deactivated  and 

drive  force  F,  =  2000  lb.  is  activated.  Restart  inte- 
b  f 

gration. 

(4)  tj  (Figure  6.24):  The  barrel  B  breaks  contact  with  the 
front  buffer  and  =  0  lb  .  Restart  integration. 

(5)  t^  (Figure  6.21):  The  barrel  B  impacts  and  captures  the 
sear  SR  which  was  locked  to  the  receiver.  The  rear 
buffer  =  12100  lb  acts  against  the  sear.  Release 
constraint  between  sear  and  receiver,  perform  momentum 
balance  with  coefficient  of  restitution  e  =  0,  activate 
constraint  between  barrel  B  and  sear  SR,  and  restart 
integration  with  new  velocities. 

(6)  tj.  (Figure  6.26):  The  barrel  B  and  sear  SR  come  to  rest. 

The  barrel  drive  force  F^  is  deactivated,  the  drive  force 

F.  is  activated,  and  the  rear  buffer  force  B  is  deacti- 
f  r 

vated.  Restart  integration. 

(7)  tg  (Figure  6.2l):  If  automatic  fire  is  to  terminate,  the 
barrel  B  and  sear  SR  return  to  the  initial  sear  position. 
The  sear  impacts  the  receiver,  and  the  sear  and  barrel 
are  captured  by  the  receiver.  Perform  momentum  balance 
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with  coefficient  of  restitution  e  =  0  and  activate  the 
constraint  between  sear  and  receiver.  The  cycle  is 
complete  with  sear  and  barrel  locked  to  receiver. 

(7')  tg  (Figure  6.2l):  If  automatic  fire  is  to  continue,  the 
barrel  B  and  sear  SR  return  to  the  initial  sear  position. 
The  sear  impacts  the  receiver  and  is  captured  by  the 
receiver,  while  the  barrel  is  released  frcm  the  sear. 
Release  the  constraint  between  sear  and  barrel,  perform 
momentum  balance  with  coefficient  of  restitution  e  =  0, 
activate  the  constraint  between  sear  SR  and  receiver,  and 
restart  integration  with  new  velocities.  The  cycle  is 
complete  and  the  barrel  is  in  the  runout  configuration 
for  another  round. 

Logical  times  t^  to  t^  depend  upon  the  state  of  the  system; 
the  relative  horizontal  displacements  and  relative  velocities  between 
bodies  of  the  system.  Since  the  horizontal  position,  velocity,  and 
acceleration  of  body  centers -of -mass  are  state  variables,  logical 
times  are  expressed  as  functions  of  these  variables. 

The  logical  events  are  defined  as  follows: 


(1) 

V 

X2 

-  34.26  =  f 

(2) 

V 

X2 

-  36.75  =  i 

(3) 

V 

X2 

-  34.26  =  i 

(4) 

V 

X2 

-  Xj  -  16  = 

(5) 

y 

X2 

=  ih  =  o 

(6) 

V 

X2 

-  Xj  -  16  = 
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The  six  events  to  tg  are  thus  defined  by  the  four  logical 
1  4 

variables  £  to  £  .  In  order  to  incorporate  these  event  predictors 
into  the  numerical  integration  algorithm,  the  derivatives  of  the 
above  equations,  with  appropriate  initial  conditions,  axe  formulated 
and  integrated  along  with  the  system  equations  of  motion.  Thus, 


;1 

t  =  x2  , 

Ao) 

=  -  18.26 

cv 

•X 

II 

cu 

£2(0) 

=  -  20.75 

,3  • 

X  —  Xp  "  X-2  y 

£3{o) 

=  0 

£h(  o) 

=  40 

The  procedure  for  determining  the  complete  system  state 

precisely  at  logical  times  to  tg,  identified  by  logical  variables 
1  4 

£  to  £  ,  is  as  follows.  An  appropriate  time  step  is  determined  by 
the  numerical  integration  algorithm  based  on  the  previous  system 
state,  polynomial  predictor  order,  and  error  tolerance.  Each  logical 
variable  in  succession  is  predicted  ahead  in  time,  using  this  time 
step.  If  no  logical  variable  is  found  to  have  passed  through  zero, 
the  program  advances  the  solution  by  the  desired  time  step  and  the 
process  is  repeated.  If  one  or  more  logical  variables  have  passed 
through  zero,  the  precise  times  at  which  the  corresponding  logical 
variables  are  zero  are  calculated  by  interpolation,  using  the  poly¬ 
nomial  predictor.  A  solution  is  then  forced  at  the  earliest  logical 
time,  indicating  occurance  of  the  first  event.  Control  is  then 
returned  to  user  supplied  subroutines  so  that  actions  can  be  taken 
according  to  the  intent  of  the  active  logical  variable. 


Discontinuous  events  can  be  categorized  according  to  the  order 
of  increasing  difficulty  as  follows : 

(1)  Those  requiring  restart  of  the  integration  process  only, 
such  as  when  discontinuous  forces  act  on  or  within  the 
system.  These  forces  are  not  considered  to  be  impulsive 
in  nature,  thus  only  discontinuous  accelerations  result. 

(2)  Those  requiring  momentum  balance  due  to  impulsive  external 
loads  (impact  between  bodies  and  mass  capture  or  release 
is  excluded)  with  no  supplemental  restitution  equations 

or  constraint  equation  modification. 

(3)  Those  requiring  momentum  balance  due  to  impact  between 
bodies  and  mass  capture  or  release.  Supplemental 
coefficient  of  restitution  equations  are  appended  to  the 
momentum  balance  equations  to  achieve  the  desired 
velocity  changes.  Constraints  are  added  or  deleted,  as 
needed  to  facilitate  mass  capture  or  release. 

The  six  events  t  to  t^  fall  into  the  following  three 
1  o 

categories : 

(1)  t  ,  t.,,  t  -  These  events  define  discontinuous  forces  of 

1  0  p 

relatively  small  magnitude,  therefore  only  a  restart  of 
the  integration  procedure  is  required. 

(2)  t  -  This  event  defines  an  externally  applied  impulsive 
load  requiring  a  momentum  balance  and  restart  of  the 
integration  procedure. 
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(3)  t^,  tg  -  These  events  define  impulsive  loading,  due  to 

impact  between  bodies  of  the  system,  and  mass  capture  and 
release.  Supplemental  equations  are  required  for  momentum 
balance  and  a  restart  of  the  integration  procedure  is 
required. 

The  effects  of  the  various  events  at  logical  times  t  to 

1  o 

on  the  position,  velocity,  and  acceleration  of  the  barrel  are  shown 
in  Figures  6.27  to  6.29,  respectively. 

The  logical  events  are  clearly  marked  by  discontinuities  in 
acceleration  of  the  gun  tube  as  shown  in  Figure  6.29.  At  t^,  an 
impulse  of  -  880  lb^-sec  is  applied  to  the  barrel,  resulting  in  a 
discontinuity  in  velocity  and  a  change  in  barrel  direction.  The 
applied  loads  are  changed  at  this  instant,  resulting  in  a  net  change 
in  acceleration. 

Impact  between  barrel  and  sear  occurs  at  t^,  resulting  in  a 
second  discontinuity  in  barrel  velocity.  Buffer  action  on  the 
combined  sear  and  barrel  mass  results  in  an  acceleration  of 
4500  in/sec  in  the  time  interval  t^  to  t,_.  The  dynamic  response 
curves  of  Figures  6.28  and  6.29  clearly  indicate  the  intermittent 
nature  of  weapon  mechanism  motion.  It  is  emphasized  that  one 
standard  set  of  equations  of  motion  and  one  standard  set  of  momentum 
balance  equations  are  generated  by  the  program.  At  the  various 
logical  times,  one  simply  appends  or  removes  constraints,  restitution 
equations,  forces,  impulses,  etc.,  as  directed  by  programmed  control 
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logic,  determined  by  the  active  logical  variables.  Thus,  an  entire 
simulation  can  be  completed  in  a  single  computer  run. 
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CHAPTER  7 

CONCLUSIONS  AND  RECOMMENDATIONS 

7.1  Conclusions 

In  this  report,  two  methods  are  presented  to  improve  efficiency 
of  a  general  purpose  dynamic  analysis  and  design  system  computer 
program.  The  first  method  employs  Gaussian  elimination  in  the  con¬ 
straint  Jacobian  matrix  to  identify  independent  and  dependent 
generalized  coordinates,  thus  defining  a  minimal  set  of  differential 
equations  whose  solution  yields  total  system  response.  The  second 
method  employs  pieced  interval  analysis  with  a  logical  events  monitor 
and  momentum  balance  capability  to  efficiently  handle  intermittent 
motion  problems. 

The  analysis  techniques  have  demonstrated  an  order  of  magni¬ 
tude  improvement  in  program  performance .  Since  these  methods  require 
no  change  in  equation  formulation,  generality  of  the  program  is 
maintained. 

The  methods  developed  are  effectively  applied  to  kinematic, 
kinetostatic ,  and  dynamic  analysis  of  constrained  planar  and  three- 
dimensional  systems.  In  addition,  current  investigations  have 
demonstrated  that  performance  of  an  analysis  program  that  includes 
body  flexibility  is  significantly  improved  when  a  minimum  number  of 
equations  of  motion  are  solved. 
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7 . 2  Recommendations 

The  method  of  coordinate  partitioning  has  potential  in 
other  areas  where  coordinates  are  related  by  algebraic  constraint 
equations.  For  example,  Heltne  [  59  J  makes  the  following  statements 
in  his  thesis  when  discussing  techniques  for  solving  large-scale 
nonlinear  programming  problems: 

(1)  "Luenberger  [  60  J  in  his  proof  of  the  convergence  rate 
of  the  GRG  algorithm  [  59  J  shows  this  convergence  rate  to  be 
coordinate  dependent.  That  is,  the  convergence  rate  is  a  function 
of  which  variables  are  selected  to  be  basic  (dependent)  and  which 
nonbasic  ( independent ) . " 

(2)  "A  second  theoretical  problem  is  the  use  of  (matrix) 
structural  considerations  for  selecting  the  basis.  For  example,  the 
use  of  a  relation  such  as  maximizing  the  trace  or  the  product  of  the 
diagonal  elements  of  a  matrix." 

He  further  states  that,  "a  nonzero  diagonal  does  not  imply  a 
nonsingular  matrix,  but  some  relation  based  upon  the  elements  of  the 
Jacobian  matrix  may  be  a  feasible  alternative." 

The  method  of  Gaussian  elimination  with  full  row  and  column 
pivoting  for  coordinate  partitioning  accomplishes  ail  of  the  above 
requirements.  Furthermore,  as  discussed  in  Chapter  3,  the  resulting 
matrix  will  generally  have  the  largest  determinant  of  all  possible 
submatrices,  it  will  generally  be  well  conditioned,  and  it  has  the 
largest  elements  on  the  diagonal. 


The  coordinate  partitioning  method  has  pitfalls  -when  the 
constraint  Jacobian  matrix  is  sparse.  To  maintain  program  efficiency 
sparse  matrix  algorithms  employ  a  partial  pivoting  strategy  that 
generally  does  not  select  the  largest  pivotal  elements.  Thus, 
coordinate  partitioning  by  this  method  will  not  be  optimum.  This 
may  cause  problems  when  attempting  to  maintain  error  control  on 
dependent  variables. 

The  approach  taken  in  this  study  is  to  periodically  construct 
a  full  matrix  version  of  the  constraint  Jacobian  and  factor  using 
the  algorithm  of  Appendix  A.  This  step,  requiring  O(r^)  operations 
and  m  x  3n  memory  locations,  is  undesirable  for  large  sparse  systems. 
Following  this,  the  sparse  matrix  is  then  augmented  to  force  zero 
changes  in  independent  variables  at  each  iteration  step.  When 
systems  experience  large  angular  displacements  the  independent 
variable  set  changes  frequently  and  correspondingly  the  sparse  matrix 
structure  changes.  Costly  symbolic  sparse  matrix  refactorizationc 
are  then  required  [  6l  ] . 

Heltne  [  59  ]  has  recently  developed  sparse  matrix  manipulation 
algorithms  employing  block  lower  triangular  factorization.  This 
technique  is  also  employed  by  the  sparse  matrix  code  in  DADS. 

However,  he  has  further  developed  an  algorithm  for  ordering 
matrices  after  a  rank-one  update.  That  is,  when  the  matrix  has  been 
put  into  block  lower  triangular  form  and  the  independent  and  dependent 
sets  change,  generally  a  new  symbolic  factorization  can  be  obtained 
at  a  fraction  of  the  cost  of  a  complete  matrix  refactorization.  He 
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has  demonstrated  an  order  of  magnitude  reduction  in  computer  cost. 
When  the  variable  sets  change  often,  which  is  the  case  when  bodies 
experience  large  angular  displacements,  the  savings  will  be  signif¬ 
icant  . 

Additional  research  required  in  this  and  related  areas  is: 

(1)  Investigate  the  most  efficient  methods  of  Gaussian 
elimination  employing  sparse  matrix  techniques  solely 
for  the  purpose  of  identifying  optimum  coordinate 
partitioning,  thus  eliminating  full  matrix  factorization. 

(2)  Implement  Heltne's  matrix  ordering  rank -one  update 
algorithm  into  the  sparse  matrix  code. 

(3)  Investigate  methods  for  determining  when  Gaussian 
elimination  is  required  for  the  purpose  of  redetermining 
the  partitioning  of  variables  into  dependent  and  inde¬ 
pendent  sets. 

(4)  Investigate  methods  of  analyzing  coordinate  velocities 
and  accelerations  as  an  alternate  procedure  for  deter¬ 
mining  the  need  for  applying  Gaussian  elimination  in  3. 

Or  as  a  complete  replacement  for  3.  One  such  technique 
might  be  to  compare  differences  between  extrapolated 
dependent  variables  (provided  by  the  integration 
algorithm)  and  the  corrected  dependent  variables  (pro¬ 
vided  by  Newton  iteration  to  satisfy  constraints).  A 
large  discrepancy  between  a  predicted  and  a  corrected 
dependent  variable  is  an  indication  that  it  should  become 
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independent.  Correspondingly  an  independent  variable 
having  a  reduced  or  minimum  velocity  and  an  acceleration 
indicating  that  it  may  remain  small  for  some  time  should 
then  be  made  dependent.  The  reverse  process  of  looking 
at  maximum  dependent  velocity  and  corresponding  accel¬ 
eration  could  also  be  employed  to  identify  a  variable 
for  the  independent  set. 

(5)  Investigate  methods  for  relating  constraint  closure 
tolerance  to  numerical  integration  error  tolerance. 

Tight  constraint  tolerance  results  in  excessive  Newton 
iterations  with  no  gain  in  accuracy.  Loose  constraint 
tolerance  results  in  poor  dependent  variable  prediction 
which  again  results  in  additional  corrector  iterations 
or  excessive  dependent  variable  error.  Methods  that 
dynamically  adjust  constraint  tolerance  to  minimize  the 
number  of  corrector  iterations  should  be  investigated. 
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APPENDIX  A 

RANK  DETERMINATION  AND  DECOMPOSITION 
OF  SINGULAR  MATRICES 

A . 1  Introduction 

The  technique  and  basic  structure  of  the  subprogram  described 
here  are  patterned  after  the  subroutine  MFGR  of  the  Scientific 
Subroutine  Package,  an  IBM  application  program  [22].  Consider  a 
system  of  equations  of  the  form 

A/  \  *  x  /  n \  =  b /  \  (A.l) 

(m  x  n)  (n  x  1)  (m  x  l) 

where  in  general,  the  matrix  A  may  lack  full  row  and/or  column  rank. 
The  task  is  to  determine  if  one  or  more  solutions  to  Equation  A.l 
exist,  and  to  find  at  least  one  efficiently  and  accurately. 

The  following  calculations  can  be  performed  on  a  general  rec¬ 
tangular  matrix  which  allows  one  to  obtain  solutions  to  Equation  A.l 
if  such  exist: 

(l)  Determine  matrix  rank  and  linearly  independent  rows  and 
columns . 


(2)  Express  a  submatrix  of  maximal  rank,  as  a  product  of 
triangular  factors. 
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(3)  Express  nonbasic  rows  in  terms  of  basic  rows. 

(U)  Express  basic  variables  in  terms  of  free  variables. 

A.  2  Theoretical  Background 

Basic  variables,  taken  here  as  dependent  variables,  will  be 
denoted  by  the  vector  u  and  free  variables,  taken  as  independent 
variables,  will  be  denoted  by  the  vector  v. 

Calculation  (l)  is  most  critical.  Matrix  rank  is  determined 
using  the  standard  Gaussian  elimination  technique,  with  complete 
pivoting.  Therefore,  the  rows  and  columns  of  the  m  by  n  matrix  A 
are  generally  interchanged  at  each  elimination  step  to  bring  the 
largest  element  to  the  pivot  position.  The  interchange  information 
is  recorded  in  two  integer  permutation  vectors  LROW  and  ICOL.:  The 
i  row  (column)  of  the  interchanged  matrix  corresponds  to  the  IROW 
(l)th  row  (ICOL  (l)th  column)  in  the  original  matrix  A.  Initially, 
IROW(J)  =  J  and  ICOL(j)  =  J. 

The  notation  A^  is  used  for  the  interchanged  matrix  implied 
fch 

at  the  i  elimination  step.  Superscripts  do  not  mean  powers.  They 
indicate  the  current  elimination  step  at  which  the  result  is  obtained. 

A. 2. a  First  Elimination  Step 

Let  a  be  the  absolutely  greatest  element  of  matrix  A,  which 
Jk 

is  found  first  in  a  columnwise  scan.  The  internal  tolerance  TOL  is 

set  equal  to  |EPS  •  a..  J.  The  parameter  EPS  is  a  specified  error 

J  ^ 

level  whose  magnitude  should  be  of  the  order  of  the  existing  computer 


round  off  error  level. 
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If  la.,  I  <  TOL,  further  calculation  is  bypassed.  Otherwise 
1  jk1  = 

rows  1  and  j  and  columns  1  and  k  of  matrix  A  are  interchanged,  giving 
a\  The  same  interchanges  must  be  applied  to  IROW,  interpreted  as  a 
column  vector,  and  to  ICOL,  interpreted  as  a  row  vector. 

A"*-  is  uniquely  expressed  as  the  product  I/*-  *  D1  •  IT*-  by 
imposing  the  following  conditions: 

(i)  L'1'  is  the  m  by  m  identity  matrix,  except  for  the  first 

column.  The  first  diagonal  element  has  a  value  of  one. 

(Il)  D1  is  an  m  by  n  matrix  with  first  diagonal  element  equal 
to  one,  while  all  remaining  elements  of  the  first  row  and 
column  are  equal  to  zero. 

(ill)  U1  is  the  n  by  n  identity  matrix,  except  for  the  first 


row. 


More  explicitly 
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(A. 2) 


Capital  le  ,irs  are  used,  as  notation  for  matrices  in  this 
appendix.  Small  letters  represent  scaler  or  vector  quantities.  The 
symbol  0  means  zero  matrices  and  I  means  identity  matrices,  where 
dimensions  are  implied  by  compatibility.  The  terms  in  Equation  A. 2 
are  further  described  as  follows: 

The  elements  of  the  first  column  of  L1  are 


l 


1 

11 


=  1, 


^4l/all’•••, 


i 

ii 


The  elements  of  submatrix  of  D~  are 


,1  1  a1 

d..  =  a,,  -  £. . 
ik  ik  il 


Ulk  =  aik 


1  .  il  Ik 


11 


i  =  2, 
k  =  2, 


. .  m 


The  elements  of  the  first  row  of 


U1 


1111  11 
U11  ail'  U12  “  a!2’  *•'’  Uln  ain 


are 
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Note  that  it  in  possible  to  record  all  nontrivial  entries  ef  1,  ,  I), 
and  U1  compactly  in  the  storage  locations  occupied  by  the  original  A, 
storing  only: 


4  ^2 


4  4 


A.2.b.  Second  Elimination  Step 


Let  d^k,  ^  =  g  be  the  absolutely  largest  element  of  D^.  If 


i&jjJ  £  TOL,  D^2  is  interpreted  as  being  the  zero  matrix. 

If  Dg2  is  not  zero  in  the  above  sense,  it  may  be  decomposed 
analogously.  In  the  compact  scheme  of  above,  rows  2  and  j ,  and  columns 
2  and  k  are  interchanged,  obtaining 


1  tt2 
U11  U12 


2  2 
L21  D22 


The  same  interchanges  must  be  performed  with  IROW  (interpreted  as  a 
column  vector),  with  ICOL  (interpreted  as  a  row  vector),  and  with  A1 
giving 
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■where 


A.2.C.  Final  Result  of  Elimination  Process 

2 

At  the  next  elimination  step  D is  factorized,  and  so  on. 
x* 

Now  assume  that  finally  D^+1  r+1  equals  zero  in  the  sense  that  all 

its  elements  are  absolutely  less  than  TOL.  This  means  that  A  has 

r  r  r  x 

the  rank  r  and  the  end  result  is  the  factorization  A  =  L  •  D  *  IT  ; 
that  is. 
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U  = 


11 


u 


u. 


12 

2 

22 


0 


1 r 


.rr 


UR  = 


The  matrix  L  is  of  dimension  r  by  r  and  unit  lower  triangular  and  U 
is  of  dimension  r  by  r  and  upper  triangular.  The  matrix  LR  is  of 
dimension  m  -  r  by  r;  if  the  given  matrix  A  is  row  regular  (that  is, 
m  =  r)  LR  is  absent  in  the  final  factorization.  The  matrix  UR  is 
of  dimension  r  by  n  -  r;  if  the  given  matrix  A  is  column  regular 
(that  is,  n  =  r)  UR  is  absent  in  the  final  factorization. 

A.2.d.  Further  Calculations  Performed 
The  problem  of  matrix  factorization  arises  in  connection  with 
the  solution  of  systems  of  equations  Ax  =  b.  Three  different  cases 
must  be  distinguished,  as  follows: 

(l)  r  =  m  =  n; 

A  is  nonsingular,  so  Ax  =  b  has  a  uniquely  determined 


solution. 


•** 
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(2)  r  <  aj 

A  is  not  row  regular,  so  solutions  of  Ax  =  b  exist  only 
if  linear  combinations  among  the  rows  of  A  are  also  valid 
among  the  rows  of  b 
(5)  r  <  n; 

A  is  not  column  regular,  so  Ax  =  0  has  nontrivial 
solutions . 

The  cases  (2)  and  (3)  may  occur  combined. 

The  solution,  if  it  exists,  is  uniquely  determined  if  r  =  n. 

Otherwise,  it  contains  n  -  r  free  parameters.  Often  one  requires  the 

linear  combinations  among  the  rows  of  a  given  matrix  A  and  the  linear 

forms  expressing  basic  variables  in  terms  of  free  ones.  Therefore, 

instead  of  IK  and  UR,  matrices  C  and  H  are  returned,  containing 

linear  combinations  and  homogeneous  solutions  respectively. 

Observe  that  the  above  calculated  factorization  belongs  to  the 
r  r  r  r 

interchanged  matrix  A  .  Therefore,  Ax  =  b  is  dealt  with  instead 

v  TOOT  ( 

of  Ax  =  b,  where  [*  J  is  obtained  from  [*],  using  the  L  '  ..J 

br  th  b  IROW(j)th 

element  of  [*]  as  the  l  element  of 

J 

r  _  . 

x  I  =  1,  . . . ,  n 
r  ’ 

b  J  ®  1,  . . . ,  m 

r  r  u  b^ 

Let  x  ,  b  be  partitioned  into  [  J  and  [^2],  Then, 

r{y  •  to'  w]  •  0  - 
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or  explicity, 


L*U*u+L*UR  •v  =  'b1 
I£*U’u+LR’UR*v=b2 
Since  L  and  U  are  nonsingular,  this  implies  that 

u  =  U-1  •  L_1  •  b1  -  U*1  •  UR  •  v 

-112 
LR  •  L  •  b  =  b 

For  convenience,  LR  is  replaced  by  LR  •  L  1  =  C  and  UR  is 

replaced  by  -  U  1  •  UR  =  H,  while  L  and  U  remain  untouched.  Consis- 

2  1 

tency  requires  that  b  =  C  •  b  ,  and  homogeneous  solutions  are  given 
by  u  =  H  •  v.  In  case  of  a  consistent  system  of  equations  Ar  •  xr 
=  br,  the  general  solution  is  xr  =  [^],  with  u=U1,L^*b+H 
•  v,  while  the  values  of  the  free  variables  contained  in  v  may  be 
chosen  arbitrarily. 


A.2.e.  Format  of  Results 

The  subprogram  returns  matrices  L,  U,  C,  H,  and  D  in  the 
storage  area  originally  occupied  by  the  input  matrix  A,  in  the  com¬ 
pact  scheme  shown  in  Figure  A.l. 


-3BV»JFry 
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